function [XM_OUT,XOW_OUT,BAI3,BAI4,Ts,Tsdry,T,CumFrac,IncFrac,Fc,sXOW_OUT, mgratio,NaKratio,test, initialBulkComp, CL_OUT,CS_OUT,CR_OUT, CL_BAI3,CS_minerals_OUT,Do_out,...
    P_out,Dn_Test,ln, WnS, Ln, Qn, Lnr, Lnexp, Wnr,Gar2SpTransition_index, Sp2PlagTransition_index,Gar2SpTransition_pressure, Sp2PlagTransition_pressure,...
    Lnexp_IT,DUTh_out, KDs,dFdP,helpDEBUG] = ...
    PETROGEN2020(Tb,Pb,dFdP,ZO,rFrac_crit,dFdPwaterchange,initialwatercontent, Co,H2O_index,Dmatrix,U_index, Th_index,IT_extended, initialTransition,forceSpPlag)
%==========================================================================
%  PETROGEN2019 Model for MORB melting in the spinel,plagioclase, and garnet fields



addpath(genpath(pwd));
% Inputs:

% Tb- Real temperature (i.e., adiabatic gradient must be added) (C)
% Pb- Pressure (kbars)
% Note: Tb & Pb must be vectors of the same length (both values
% should be either constant or monotonically decreasing as if
% ascending along an adiabat).
% dFdP - Melt productivity (melt fraction/kbar).  Can either be
%           specified either as a single value or vector of length Tb & Pb.
% ZO - Mantle bulk composition in [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;
% rFrac_crit - retained melt fraction during melting where 0 = pure fractional
%               melting and 1 = pure batch melting
%dFdPwaterchange - factor that dFdP is divided by when melting occurs due
%                   to the presence of water (10 means that dFdP is reduced
%                   by a factor of 10)
%initialwatercontent - initial water contnent above nominally anhydrous
%conditions in ppm
%Co_Designation  - Initial trace element source composition
%                     1 = Depelted MORB Mantle from Workman and Hart (2005)
%                     2 = Primitive Upper Mantle from Hart and Zindler (1986)
%                     3 = Initial Garnet Signature (calcuated using PETROGEN)
%Dmatrix  = matrix of D's that match the elements in Co_Designation with columns of minerals
%                [Aug, Opx, Oliv, Plag, Sp]
%  IT_extended = 2D temperature field
% initialTransition - parameter neccessary to calculate initial phase
% transition pressures, odn't worry about changing this




%  Outputs
% XM_OUT - Mantle mode [Aug, Opx, Oliv, Plag, Sp]
% XOW_OUT - Mantle composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% BAI3 - All incremental melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% BAI4 - Extracted incremental melts [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% Ts - Solidus temperature (C)
% Tsdry - Dry solidus temperature (C)
% T - Temperature corrected for latent heat of melting based on the
%       amount of melt extracted at each step (C)
% CumFrac, IncFrac - cumulative and incrmental melt fraction
% Fc  - max(CumFrac)
% rXOW - Residual mantle composition (residue = retained melt + solid) [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] (wt%)
% MgRatio -  [Mg# residue, Mg# melt] ;
% NaKratio - NaK# melt
% test - [Pb Ts Tb (1-XMGliq) WNA WT WK WNACA Sp2PlagTransitionPressure Gar2SpTransitionPressure QCA];
% initialBulkComp - the mantle bulk composition used
% CL_OUT - trace element composition of extracted incremental melts
% CS_OUT - trace element composition of the solid
% CR_OUT - trace element composition of the residue
%           (residue = retained melt + solid)
% CL_BAI3 - trace element composition of all incremental melts
% CS_minerals_OUT - trace element composition minerals in solid residue, right now it is just CPX
% Do_out - the bulk modal-melting D at each step
% P_out - the non-modal melting P at each step
% Dn_Test - the non-modal melting D at each step
% ln - mass of new liquid generated at each step
% WnS - mass of the solid in the residue at each step
% Ln - mass of the liquid at each step
% Qn - fraction of the melt remains in the residue at each step
% Lnr - mass of melt that remains in the residue  at each step
% Lnexp - mass of melt expelled at each step
% Wnr - mass of residue at each step
% Gar2SpTransition_index - index of the garnet to spinel phase transition
%                           (i.e, phase transition pressure = Pb(Gar2SpTransition_index)
% Sp2PlagTransition_index - index of the spinel to plagioclase phase transition
%                           (i.e, phase transition pressure = Pb(Sp2PlagTransition_index)

% Lnexp_IT -  mass of melt expelled in 2D
% DUTh_out -  the bulk DU/DTh at each step (which is pressure dependent)


%
%  Example of batch melting with HZ dep 1 mantle composition:
%     [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
%            BG15([1350,1340],[12,11],dFdP,2,2,0.01);

%  Example of incremental batch melting with 90% melt extraction at each
% step using the HZ dep 1 mantle composition:
%  [XM,XOW,BAI,BAI3,Ts,T,CumFrac,IncFrac,Fc] = ...
% BG15([1350,1340],[12,11],dFdP,2,2,0.9);
%
%  PETROGEN Version 1.0
%  Stephanie M. Brown, Mark Behn & T.L. Grove [JGR, Submitted, January 2015]
%==========================================================================

%% Initialization Sequence
% Define Oxide Weights

%[SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,MnO,CaO,K2O,Na2O,Fe2O3,P2O5];
FK = [60.0840,79.8988,101.9612,151.9902,71.8464,40.3044,70.9374,56.0794,...
    94.2036,61.9790,159.6920,141.945];

%Initial Guesses for phase transition:
Tb_GS = Tb(1);
Pb_GS= Pb(1);
dFdP_GS = rFrac_crit;


if ~exist('initialTransition','var')
    [XM_OUT,XOW_OUT,BAI3,BAI4,Ts,Tsdry,T,CumFrac,IncFrac,Fc,sXOW_OUT, mgratio,NaKratio,test]=...
        PETROGEN2019(Tb_GS,Pb_GS,dFdP_GS,ZO,rFrac_crit,1,initialwatercontent, Co, H2O_index,Dmatrix,U_index, Th_index,IT_extended,'initialTransition');
    Gar2SpTransition = test(10);
    Sp2PlagTransition = 9;
else
    Gar2SpTransition = 24;
    Sp2PlagTransition = 9;
end




Sp2PlagTransition_forced = 9;



% Find the appropriate phase transition indicies
if Gar2SpTransition > min(Pb) && Gar2SpTransition < max(Pb)
    tmp = abs(Pb-Gar2SpTransition);
    [idx Gar2SpTransition_index] = min(tmp); %index of closest value
    Gar2SpTransition_pressure = Pb(Gar2SpTransition_index);
else
    Gar2SpTransition_pressure = Gar2SpTransition;
end
%will be updated later:
if Sp2PlagTransition > min(Pb) && Sp2PlagTransition < max(Pb)
    tmp = abs(Pb-Sp2PlagTransition);
    [idx Sp2PlagTransition_index] = min(tmp); %index of closest value
    Sp2PlagTransition_pressure = Pb(Sp2PlagTransition_index);
else
    Sp2PlagTransition_pressure = Sp2PlagTransition;
end

if isnan(Gar2SpTransition_pressure) == 1
    Gar2SpTransition_pressure=24;
end



%Melting reaction [Aug, Opx, Oliv, Plag, Sp]
YMP = [.31,0.16,-0.06,0.59,0,0]; %Plag
YMS = [0.97,0.25,-0.30,0.00,0.08,0]; % Spin (init);  P-dependent spinel reaction is defined below
YMG = [0.81,-0.5,0.08,0,0,0.60]; % Garnet
YMG = YMG./sum(YMG);

%Normalize initial mantle bulk composition
ZO = ZO./sum(ZO);
% Calculate initial modes, needs to know the Mg# so specifify a KD
KD = 0.33;
[XMG,XMS,XMP] = calculateInitialModes(ZO,Pb(1),KD,Gar2SpTransition_pressure,Sp2PlagTransition_pressure);

% Choose Mantle Mode and Composition to be used.
XO = ZO;
XOW = XO;

% Normalize (again)
XO=100*XO./sum(XO);
XOW = 100*XOW./sum(XOW);
XMS = XMS./sum(XMS);
XMP = XMP./sum(XMP);
XMG = XMG./sum(XMG);

%Define Major Element Partition Coefficients (Note: some will be updated based on
%   pressure below).

% Columns are Ti, Cr, K, Na
% Rows are cpx, opx, olv, plagioclase/spinel/garnet
DMP = ...
    [0.42	4.3       0.001	    0.15
    0.19	5         0.001	    0.04
    0.03	1.45     0.001	   0.001
    0.001   0.001	0.271	  1];

DMS = ...
    [0.22,  4.3,     0.001,   0.25;
    0.14,   5.00,   0.001,   0.06;
    0.02,   1.45,   0.001,   0.001;
    0.14,   75.00,  0.001,   0.001];

DMG = ...
    [0.22,  5.5,     0.001,  0.25;
    0.14,   4.00,   0.001,  0.06;
    0.02,   1.45,   0.001,  0.001;
    0.30,   5.00,   0.001,  0.01];





Co(H2O_index) = initialwatercontent;
Co = Co;
% Additional parameters
LH=4e5;             %Latent Heat
Cp=1e3;             %Specific heat

% % Hess 1992
LH=195;
Cp = 0.3;

% LH=400;
% Cp =1;

% Initialize dFdP array
if length(dFdP)==1
    dFdP = ones(1,length(Pb))*dFdP;
    dFdPtype = 'constant';
else
    dFdPtype = 'variable';
end


% Initialize melt fraction
MeltFrac = zeros(1,length(Pb));  %Array for melt fractions
for i = 1:length(Pb)
    if length(Pb) == 1
        MeltFrac(i) = dFdP;
    else if size(unique(Pb),2)==1
            MeltFrac(i) = dFdP(i);
        else
            if i == 1
                MeltFrac(i) = -(Pb(i+1)-Pb(i)).*dFdP(i);
            else
                if Pb(i)-Pb(i-1)==0
                    fixPbx(i)=i;
                    MeltFrac(i) = dFdP(i);
                else
                    MeltFrac(i) = -(Pb(i)-Pb(i-1)).*dFdP(i);
                end
            end
        end
    end
end




if exist('fixPbx')==1
    MeltFrac(find(fixPbx,1)-1) = dFdP(find(fixPbx,1)-1);
    dFdPtype = 'constant';
end

% Initialize BAI, BIW, and BAI3 arrays that will hold melt compositions
BIW = zeros(1,9);   %Incremental Melt (wt%)
%rBIW = zeros(1,9);   %Retained Melt (wt%)
BAI3 = zeros(length(Pb),9);   %Incremental Melt (cation units) %XX
BAI4 = zeros(length(Pb),9);   %Extracted Incremental Melt (cation units)
CL_BAI3= zeros(length(Pb),size(Co,1));

% CS_OUT= zeros(length(Pb),size(Co,1));
% CR_OUT= zeros(length(Pb),size(Co,1));
testXMG =XMG(1:2);
Cabulk_mode = 1;
%%
% Calculate initial pyroxene modes in spinel field based on P-dependent
%   cpx-opx join
MgNumBulk = ((XO(:,6)/40.311))./(((XO(:,6)./40.311)+((XO(:,5)./71.846))));
MgLiqEstimate = 1/(1+(1/KD)*((1-MgNumBulk)/MgNumBulk));

PCA = calculatePCA(Pb(1), MgLiqEstimate, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
CaOOPX = 2;


%Spinel Field
Pyx = XMS(1) + XMS(2);
Cabulkfactor = 1./Pyx;
XMS(1) = Pyx*(Cabulkfactor*XO(7) - CaOOPX)./(PCA - CaOOPX);
XMS(2) = Pyx-XMS(1);
if XMS(1) <= 0.; XMS(1) = 0.0;  XMS(2) = Pyx; end
if XMS(2) <= 0.; XMS(2) = 0.0;  XMS(1) = Pyx;end
XMS = XMS./sum(XMS);



% Calculate initial pyroxene modes in garnet field based on P-dependent
%   cpx-opx join
clear Pyx
Pyx = XMG(1) + XMG(2);
Cabulkfactor = 1./Pyx;
PPX = Cabulkfactor*(XOW(7) - 4.96*XMG(6));
XMG(1) = Pyx*(PPX - CaOOPX)./(PCA - CaOOPX);
XMG(2) = Pyx-XMG(1);
if XMG(1) <= 0.; XMG(1) = 0.0; XMG(2) = Pyx;  end
if XMG(2) <= 0.; XMG(2) = 0.0; XMG(1) = Pyx; end
XMG = XMG./sum(XMG);

%%
% Define mantle mode arrays
XOW_OUT = zeros(length(Pb),length(XOW));
XM_OUT = zeros(length(Pb),length(XMS));

XOW_OUT(1,:) = XOW;
T = zeros(1,length(Pb));
Ts = zeros(1,length(Pb));

sXOW = XOW;
sXOW_OUT(1,:) = sXOW;
initialBulkComp = XOW;

%% Iteration Sequence
%Array for corrected temperature
%Array for solidi
iiii= 1;
Frac = 0 ;     % Melt fraction is set to 0 at first step
rFrac = 0;
CumFrac = zeros(1,length(Pb));
cumFrac = 0;
Lnexp_IT = zeros(size(IT_extended));



IncFrac = zeros(1,length(Pb));
% mass balance inital parameters
Wo = 1;
Lnr(1) = 0;
WnS(1) = Wo;
Wnr(1) = Wo;
CR = Co;
CS = Co;

IDFirstMelt=0 ;
melton = 0;         % Flag to determine whether melting is turned on.
if length(Pb) <= 2  % Note that usually melting is turned off to start,
    melton = 1;     %   but for special cases with a single melt step
end               %   it is turned on for the initial step.




% Calculate melting looping through all pressures (P)
for ii = 1:length(Pb)
    
    Frac = MeltFrac(ii);
    
    % Calculate P-dependent melting reaction [Aug, Opx, Oliv, Plag, Sp, Garnet]
    YMS(1) = 0.0016*Pb(ii)^2 + 0.0006*Pb(ii) + 0.5;
    YMS(2) = -0.0022*Pb(ii)^2 - 0.0223*Pb(ii) + 0.9365;
    YMS(3) = 0.0015*Pb(ii)^2 - 0.008*Pb(ii) - 0.3414;
    YMS(4) = 0;
    YMS(5) = 0.0024*Pb(ii) + 0.0738;
    YMS(6) = 0;
    YMS = YMS./sum(YMS);
    
    if XMG(2) > 0
        YMG = [0.81,-0.5,0.08,0,0,0.60]; % Garnet
        YMG = YMG/sum(YMG);
    else
        YMG(2) = 0;
        YMG = YMG/sum(YMG);
    end
    
    
    
    
    % Estimage Mg# of melt from Sp Lherz as function of FeO, MgO, F, Phase
    % Proportions and mineral-melt Fe-Mg Kds
    
    %Re-order elements in original KG92 formulation for BG15_Mgnum routine:
    %  [K,Na,Ca,Fe,Mg,Ti,Si,Al]
    
    % [1-SiO2,2-TiO2,3-Al2O3,4-Cr2O3,5-FeO,6-MgO,7-CaO,8-K2O,9-Na2O]
    XOW = 100.*XOW./sum(XOW);
    
    
    %update DU and DTh for CPX based on pressure using P regression on
    %Table 4 in Landwehr2001
    xxxx=Pb(ii)/10;
    if xxxx > 3
        xxxx = 3; %keep constant above 3 gpa
    end
    
    p1 = -0.0055;
    p2 = 0.0275;
    DTh_changeCPX_P = p1*xxxx + p2;
    
    p1 = 0.2;
    p2 = 0.73;
    % p1 = 0.130;
    % p2 = 0.895;
    %   p1 = 0.12983;
    %   p2 = 0.89477;
    DUDTh_changeCPX_P = p1*xxxx + p2;
    DU_changeCPX_P =DUDTh_changeCPX_P.*DTh_changeCPX_P;
    
    
    
    Dmatrix(U_index,1)=DU_changeCPX_P;
    Dmatrix(Th_index,1)=DTh_changeCPX_P;
    DUTh_out(ii,:) = [DU_changeCPX_P DTh_changeCPX_P DUDTh_changeCPX_P];
    
    
    
    
    
    % calculate trace element melting
    if Pb(ii) >= Gar2SpTransition_pressure
        Do = Dmatrix*XMG';
        P = Dmatrix*YMG';
    else if Pb(ii) >= Sp2PlagTransition_pressure
            Do = Dmatrix*XMS';
            P = Dmatrix*YMS';    
        else
            Do = Dmatrix*XMP';
            P = Dmatrix*YMP';
        end
    end
    
    
    if ii> 1
        Wnriiminus1   = Wnr(ii-1);
        ln_temp =  Frac.*WnS(ii-1);
        WnS_temp = WnS(ii-1) - ln_temp;
        Ln_temp = ln_temp + Lnr(ii-1);
        %Dn = (Do - P.*(Ln_temp))./(1-(Ln_temp));
        Dn = (Do - P.*(Frac))./(1-(Frac));
        CL = CR.*Wnr(ii-1)./(Dn.*WnS_temp+Ln_temp);
        
        
        
    else
        Wnriiminus1   = Wo;
        ln_temp = Frac.*Wo;
        WnS_temp = Wo - ln_temp;
        Ln_temp = ln_temp;
        %Dn = (Do - P.*(Ln_temp))./(1-(Ln_temp));
        Dn = (Do - P.*(Frac))./(1-(Frac));
        CL = Co.*Wo./(Dn.*WnS_temp+Ln_temp);
    end
    
    % Original Fe-Mg KD values from KG92
    % XOlK = 0.3;
    % XCpxK = 0.28;
    % XOpxK = 0.27;
    % XSpK = 0.42;
    % Olivine Fe-Mg Kd     (0.30)
    % Cpx Fe-Mg Kd   (0.98*XOlK = 0.294)
    % Opx Fe-Mg Kd   (0.96*XOlK = 0.288)
    % Spinel Fe-Mg Kd (1.52*XOlK = 0.456)
    
    
    
    
    if Pb(ii) >= Gar2SpTransition_pressure %Garnet field
        %Fe-Mg Kds
        if Pb(ii) > 27
            XOlK = 0.32;
            %XOlK = 0.33 + 0.02.*Pb(ii);
            
            
        else
            %XOlK = 0.33;
            XOlK= 0.003.*Pb(ii) + 0.245;
        end
        XCpxK = 1.04*XOlK;
        XOpxK = 0.93*XOlK;
        XGarK = 1.4*XOlK;
        XSpK = 1.5*XOlK;
        
        
        XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMG,YMG,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
        helpDEBUG(ii,:) = [XMGliq, Frac, Wnriiminus1, WnS_temp, XMG,YMG,XOlK,XCpxK,XOpxK,XSpK,XOW];
        
    else if Pb(ii) >= Sp2PlagTransition_pressure  %Spinel field
            
            
            
            %Fe-Mg Kds
            %XOlK = 0.33;
            
            XOlK= 0.003.*Pb(ii) + 0.245;
            
            XCpxK = 0.94*XOlK;
            XOpxK = 0.97*XOlK;
            XGarK = 1.4*XOlK;
            XSpK = 1.5*XOlK;
            
            %       ii
            %     size(Pb)
            %     Pb(ii)
            %     XOlK
            %     Tb(ii)
            %     [Frac, Wnriiminus1, WnS_temp, XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW]
            
            
            % Frac
            % Wnriiminus1
            % WnS_temp
            
            % YMS
            XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
            
            
            helpDEBUG(ii,:) = [XMGliq,  Frac, Wnriiminus1, WnS_temp, XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XOW];
            
            
        else % plagioclase field
            %Fe-Mg Kds
            %XOlK = 0.29;
            XOlK= 0.003.*Pb(ii) + 0.245;
            XCpxK = 0.94*XOlK;
            XOpxK = 0.97*XOlK;
            XGarK = 1.4*XOlK;
            XSpK = 1.5*XOlK;
            
            
            XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMP,YMP,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
            
            helpDEBUG(ii,:) = [XMGliq,  Frac, Wnriiminus1, WnS_temp, XMP,YMP,XOlK,XCpxK,XOpxK,XSpK,XOW];
        end
    end
    
    KDs(ii,:) = [XCpxK XOpxK XOlK  0 XSpK XGarK];
    
    % Calculate Mg# of olivine in residual
    XMGOL = (XMGliq/(XOlK*(1-XMGliq)))./(1+(XMGliq/(XOlK*(1-XMGliq))));
    
    % % % % % % Define P-dependent coefficients in D matrix
    % DMS(1,1) = 22.570*(Pb(ii)^(-1.53));    % Dticpx
    % DMS(1,4) = -0.440+0.557*log10(Pb(ii)); % Dnacpx
    % DMS(2,1) =  0.926*(Pb(ii)^(-0.565));   % Dtiopx
    % DMS(2,4) = -0.151+0.166*log10(Pb(ii)); % Dnaopx
    %
    % % Define P-dependent coefficients in D matrix
    % DMG(1,1) = 22.570*(Pb(ii)^(-1.53));    % Dticpx
    % DMG(1,4) = -0.440+0.557*log10(Pb(ii)); % Dnacpx
    % DMG(2,1) =  0.926*(Pb(ii)^(-0.565));   % Dtiopx
    % DMG(2,4) = -0.151+0.166*log10(Pb(ii)); % Dnaopx
    
    
    xxxx2=Pb(ii);
    if xxxx2 > 30
        xxxx2 = 30; %keep constant above 3 gpa
    end
    
    % Columns are Ti, Cr, K, Na
    % Rows are cpx, opx, olv, plagioclase/spinel/garnet
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGED ON APRIL 7 2020!!!! XX SMB S
    % WAS: Define P-dependent coefficients in D matrix
    DMS(1,1) =  0.34113 - 0.0056005.*xxxx2;    % Dticpx
    DMS(2,1) =   0.16726 - 0.0020215.*xxxx2;   % Dtiopx
    DMS(1,4) = 0.036569+0.012871.*xxxx2; % Dnacpx
    DMS(2,4) = 0.006879+0.0027618.*xxxx2; % Dnaopx
    
    %Define P-dependent coefficients in D matrix
    DMG(1,1) =  0.34113 - 0.0056005.*xxxx2;    % Dticpx
    DMG(2,1) =   0.16726 - 0.0020215.*xxxx2;   % Dtiopx
    DMG(1,4) = 0.036569+0.012871.*xxxx2; % Dnacpx
    DMG(2,4) = 0.006879+0.0027618.*xxxx2; % Dnaopx
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%     % % decrease TiO2 D XXX SBK PLAYING
%     DMP(:,1) = DMP(:,1)./2;
%     DMS(:,1) = DMS(:,1)./2;
%     DMG(:,1) = DMG(:,1)./2;


    
    %Plag field
    XCP = [XMP(1:4)]*DMP; %bulk partition coefficient in solid
    XPP = [YMP(1:4)]*DMP;  %bulk partition coefficient in melt
    DP = (XCP - XPP.*Frac)./(1-Frac);
    
    %Spinel field
    XCS = (DMS'*[XMS(1:3),XMS(5)]')'; %bulk partition coefficient in solid
    XPS = (DMS'*[YMS(1:3),YMS(5)]')'; %bulk partition coefficient in melt
    DS = (XCS - XPS.*Frac)./(1-Frac);
    
    
    
    %Garnet field
    XCG = [XMG(1:3),XMG(6)]*DMG; %bulk partition coefficient in solid
    XPG = [YMG(1:3),YMG(6)]*DMG; %bulk partition coefficient in melt
    DG = (XCG - XPG.*Frac)./(1-Frac);
    
    
    % Calculate concentration of oxide in melt (wt%)
    %Plagioclase
    XLFP(1) = XOW(2).*Wnriiminus1/(DP(1).*WnS_temp+Ln_temp); %Ti
    XLFP(2) = XOW(4).*Wnriiminus1/(DP(2).*WnS_temp+Ln_temp); %Cr
    XLFP(3) = XOW(8).*Wnriiminus1/(DP(3).*WnS_temp+Ln_temp); %K
    XLFP(4) = XOW(9).*Wnriiminus1/(DP(4).*WnS_temp+Ln_temp); %Na
    
    %Spinel
    XLFS(1) = XOW(2).*Wnriiminus1/(DS(1).*WnS_temp+Ln_temp); %Ti
    XLFS(2) = XOW(4).*Wnriiminus1/(DS(2).*WnS_temp+Ln_temp); %Cr
    XLFS(3) = XOW(8).*Wnriiminus1/(DS(3).*WnS_temp+Ln_temp); %K
    XLFS(4) = XOW(9).*Wnriiminus1/(DS(4).*WnS_temp+Ln_temp); %Na
    
    
    %Garnet
    XLFG(1) = XOW(2).*Wnriiminus1/(DG(1).*WnS_temp+Ln_temp); %Ti
    XLFG(2) = XOW(4).*Wnriiminus1/(DG(2).*WnS_temp+Ln_temp); %Cr
    XLFG(3) = XOW(8).*Wnriiminus1/(DG(3).*WnS_temp+Ln_temp); %K
    XLFG(4) = XOW(9).*Wnriiminus1/(DG(4).*WnS_temp+Ln_temp); %Na
    
    [PCA, DCA] = calculatePCA(Pb(ii), XMGliq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
    QCA = PCA/DCA;  %Ca content (XLFS)
    
    
    % a la:
    % D = Min / Melt
    % Melt = Min / D
    %CaO melt = CaO_Cpx/ D_CaO_cpx
    
    
    %         % if there is no opx then PCA should be set to the CaO in Pyx.
    %         % SMB XX: should we set DCA to anything in particular too in this case?
    %         if Pb(ii)>=Gar2SpTransition_pressure && XMG(2)== 0
    %             'here'
    %             Pb(ii)
    %             PCA
    %             PCA = (XOW(7) - 4.96*XMG(6))./(XMG(1) + XMG(2));
    %             PCA
    %             QCA = PCA/DCA;  %Ca content (XLFS)
    %         end
    
    % Calculate NaK# wt, Mg# mol, assign wt% TiO2, K2O, Na/Na+Ca wt
    
    tempNaK_Plag =[(XLFP(3)+XLFP(4))/(XLFP(3)+XLFP(4)+QCA)];
    
    if Pb(ii) >= Gar2SpTransition_pressure
        %Garnet Field
        WNA=(XLFG(3)+XLFG(4))/(XLFG(3)+XLFG(4)+QCA); %NaK# = (Na + K) / (Na + K +Ca)
        WNACA=XLFG(4)/(XLFG(4)+QCA); %Na# = Na / (Na + Ca)
        % assign wt.% TiO2 and K2O
        WT=XLFG(1);
        WK=XLFG(3);
        
    else if Pb(ii) >= Sp2PlagTransition_pressure
            %Spinel Field
            WNA=(XLFS(3)+XLFS(4))/(XLFS(3)+XLFS(4)+QCA); %NaK# = (Na + K) / (Na + K +Ca)
            WNACA=XLFS(4)/(XLFS(4)+QCA); %Na# = Na / (Na + Ca)
            % assign wt.% TiO2 and K2O
            WT=XLFS(1);
            WK=XLFS(3);
        else
            %Plag Field
            WNA=(XLFP(3)+XLFP(4))/(XLFP(3)+XLFP(4)+QCA);
            WNACA=XLFP(4)/(XLFP(4)+QCA);
            % assign wt.% TiO2 and K2O
            WT=XLFP(1);
            WK=XLFP(3);
        end
    end
    
    
    
    
    
    water(ii) = CL(H2O_index)*10^-4;
    mgratio(ii,:) = [((XOW(6)/40.311)/(XOW(6)/40.311+XOW(5)/71.846)), XMGliq] ;
    NaKratio(ii) = ((XOW(9)+XOW(8))/(XOW(9)+XOW(8)+XOW(7)))./WNA;
    
    
    % Calculate mineral components and T
    CaOAl2O3_ratio = 0.0099.*Pb(ii) + 0.48;
    DATAHERE = [1 Pb(ii) XMGliq WNA WT WK CaOAl2O3_ratio];
    
    
    if Pb(ii) >= Gar2SpTransition_pressure
        %Garnet Field
        
%         Original
        Ts(ii)=1320.279+8.468*Pb(ii)-203.208*(1-XMGliq)-131.937*WNA+32.949*WT;        
        S(1)=0.152+0.006*Pb(ii)+0.245*(1-XMGliq)-0.214*WNA+0.001*WT;  %Olivine
        S(2)=0.190+0.000*Pb(ii)+0.074*(1-XMGliq)-0.114*WNA-0.014*WT;  %CPX
        S(3)=0.727-0.008*Pb(ii)-0.343*(1-XMGliq)+0.323*WNA-0.008*WT;  %Plag
        S(4)=-0.028+0.002*Pb(ii)-0.002*(1-XMGliq)-0.397*WNA+0.018*WT; %Qtz
        
        
        
        % RevPet
%         garnet_coefficients = [...
%             1136.36052	8.73915	184.88412	-19.48245	29.00187	-23.41730	-22.48205
%             0.53105	0.00664	-0.31382	-0.18522	-0.01025	-0.02342	-0.10713
%             -0.19232	-0.00305	0.10861	0.35773	-0.00495	-0.01078	0.41337
%             0.61340	-0.00638	0.23735	0.36467	-0.01798	-0.04181	-0.24342
%             0.04833	0.00275	-0.02572	-0.54685	0.01868	0.01032	-0.06905
%             -77.53423	139.87872	46.43963	28.72057	1.52179	3.25095	18.62211];
%         
%        %Ts(ii) = DATAHERE*garnet_coefficients(1,:)';
%         for sss2=1:4
%             S(sss2)=DATAHERE*garnet_coefficients(sss2+1,:)';
%         end
      

        Tsdry(ii) = Ts(ii);
        Sdry =S;
        Ts(ii) = Ts(ii) - 26.3619.*water(ii);
        S(1)=S(1)-0.0029.*water(ii);%Olivine
        S(2)=S(2)-0.0093.*water(ii); %CPX
        S(3)=S(3)+0.0075.*water(ii);  %Plag
        S(4)=S(4)+0.0068.*water(ii); %Qtz
        
    else if Pb(ii) >= Sp2PlagTransition_pressure    %Spinel Field
            
            %   %normalized to the 7 tormey mineral components
            Ts(ii)=1212.464+11.966*Pb(ii)-96.899*(1-XMGliq)-89.350*WNA+2.367*WT; %MITCMAS
            S(1)=0.125+0.008*Pb(ii)+0.197*(1-XMGliq)-0.260*WNA-0.023*WT;  %Olivine % All Spinel experiments
            S(2)=0.181+0.000*Pb(ii)+0.071*(1-XMGliq)-0.255*WNA+0.005*WT;  %CPX
            S(3)=0.516-0.001*Pb(ii)-0.130*(1-XMGliq)+0.782*WNA-0.035*WT;  %Plag
            S(4)=0.149-0.005*Pb(ii)-0.098*(1-XMGliq)-0.379*WNA+0.025*WT; %Qtz % All Spinel experiments
            
%             spinel_coefficients =   [...
%                 1049.24929	12.71243	63.47484	-3.32516	2.65802	-12.03128	117.75713
%                 0.29432	0.00879	-0.21515	-0.20805	-0.02668	-0.00335	0.04626
%                 -0.25601	-0.00049	0.01598	0.33606	-0.00051	-0.00570	0.52520
%                 0.82612	-0.00360	0.11628	0.39210	-0.01170	-0.07164	-0.48226
%                 0.14785	-0.00478	0.08470	-0.53037	0.02331	0.01994	-0.10667
%                 -29.49925	84.81760	24.99866	24.66609	2.78588	-0.13748	-1.84832];
%             
%            % Ts(ii) = DATAHERE*spinel_coefficients(1,:)';
%             for sss2=1:4
%                 S(sss2)=DATAHERE*spinel_coefficients(sss2+1,:)';
%             end
%         
            
            Tsdry(ii) = Ts(ii);
            Sdry =S;
            Ts(ii) = Ts(ii) - 26.19.*water(ii);
            S(1)=S(1)-0.0029.*water(ii);%Olivine
            S(2)=S(2)-0.0093.*water(ii); %CPX
            S(3)=S(3)+0.0075.*water(ii);  %Plag
            S(4)=S(4)+0.0068.*water(ii); %Qtz
        else
            %Plag Field
            
            %normalized to the 7 tormey mineral components
            Ts(ii)=1215.133+10.576*Pb(ii)-72.533*(1-XMGliq)-198.174*WNA+23.590*WT;
            S(1)=0.134+0.006*Pb(ii)+0.171*(1-XMGliq)-0.249*WNA-0.012*WT;  %Olivine
            S(2)=0.230-0.003*Pb(ii)-0.057*(1-XMGliq)-0.139*WNA-0.006*WT;  %CPX
            S(3)=0.390+0.020*Pb(ii)-0.075*(1-XMGliq)+0.485*WNA-0.032*WT;  %Plag
            S(4)=0.225-0.018*Pb(ii)-0.003*(1-XMGliq)-0.229*WNA+0.009*WT; %Qtz
            
%             plagioclase_coefficients = ...
%                 [1074.38633	11.86431	65.55420	-138.22714	20.55173	5.85532	79.01883
%                 0.41221	0.00526	-0.16238	-0.36336	-0.01093	-0.00063	-0.12587
%                 -0.25604	0.00036	0.03202	0.33938	-0.00408	-0.00956	0.49666
%                 0.42923	0.01181	0.13163	0.50259	0.00140	-0.07267	-0.15741
%                 0.42247	-0.01757	0.00135	-0.48658	-0.00104	0.02298	-0.22459
%                 -43.58532	136.94746	24.54202	37.54688	2.10967	-0.77250	1.63584];   
%             % Ts(ii) = DATAHERE*plagioclase_coefficients(1,:)';
%             for sss2=1:4
%                 S(sss2)=DATAHERE*plagioclase_coefficients(sss2+1,:)';
%             end
%             
            Tsdry(ii) = Ts(ii);
            Sdry =S;
            Ts(ii) = Ts(ii) - 26.19.*water(ii);
            S(1)=S(1)-0.0029.*water(ii);%Olivine
            S(2)=S(2)-0.0093.*water(ii); %CPX
            S(3)=S(3)+0.0075.*water(ii);  %Plag
            S(4)=S(4)+0.0068.*water(ii); %Qtz
        end
    end
    
    % updates pressure of phase boundaries:
    testPboundary_Sp2Gar = (1.2  - .52*WNA+2.11*XMGliq)*10;
    testPboundary_Plag2Sp = (0.823  +1.138*WNA+0.165*XMGliq)*10;
    testPboundary_Plag2Sp_PlagSide = (0.823+1.138*tempNaK_Plag+0.165*XMGliq)*10; %This is using the Mg\# from the spinel field
    
    % if ii == 701
    %     testPboundary_Plag2Sp
    %     testPboundary_Plag2Sp_PlagSide
    %
    %     tempNaK_Plag
    %     XMGliq
    %
    %     return
    % end
    
    if testPboundary_Plag2Sp_PlagSide < testPboundary_Plag2Sp
        testPboundary_Plag2Sp = testPboundary_Plag2Sp_PlagSide;
    end
    
    
    
    if Pb(ii) >= Gar2SpTransition_pressure
        
        if all(testPboundary_Sp2Gar >= min(Pb) & testPboundary_Sp2Gar <= max(Pb))==1
            tmp = abs(Pb-testPboundary_Sp2Gar);
            [idx Gar2SpTransition_index_temp] = min(tmp); %index of closest value
            Gar2SpTransition_pressure = Pb(Gar2SpTransition_index_temp) ;
            if (Pb(ii) - Gar2SpTransition_pressure) < 0
                Gar2SpTransition_pressure = Pb(ii);
            end
        else
            Gar2SpTransition_pressure=testPboundary_Sp2Gar;
        end
        
    else if Pb(ii) >= Sp2PlagTransition_pressure    %Spinel Field
            
            if all(testPboundary_Plag2Sp >= min(Pb) & testPboundary_Plag2Sp <= max(Pb))==1
                tmp = abs(Pb-testPboundary_Plag2Sp);
                [idx Sp2PlagTransition_index_temp] = min(tmp); %index of closest value
                Sp2PlagTransition_pressure = Pb(Sp2PlagTransition_index_temp);
                if (Pb(ii) - Sp2PlagTransition_pressure) < 0
                    Sp2PlagTransition_pressure = Pb(ii);
                end
            else
                Sp2PlagTransition_pressure = testPboundary_Plag2Sp;
            end
            
            
        end
        
    end
    
    
    
    
    if strcmp(forceSpPlag,'force Sp-Plag Boundary')==1
        Sp2PlagTransition_pressure = Sp2PlagTransition_forced;
    end
    
    
    
    S = S./sum(S);
    Sdry = Sdry./sum(Sdry);
    
    
    
    test(ii,:) = [Pb(ii) Ts(ii) Tb(ii) (1-XMGliq) WNA WT WK WNACA testPboundary_Plag2Sp testPboundary_Sp2Gar QCA Sp2PlagTransition_pressure Gar2SpTransition_pressure];
    testS(ii,:) = [S sum(S)];
    testSdry(ii,:) = Sdry;
    
    
    
    % CHECK FOR MELTING
    T(ii) = Tb(ii) - ((cumFrac).*LH./Cp);
    % LH=4e5;             %Latent Heat
    % Cp=1e3;
    % cumFrac
    
    % Cp = .3-.33 cal/gm/K
    % LH=195 cal/gm
    
    
    
    % iterate and recalculate if melting occurs because of water
    
    if T(ii) < Tsdry(ii) && water(ii) > 0
        
        MeltFrac(ii) = MeltFrac(ii)./dFdPwaterchange;
        Frac = MeltFrac(ii);
        
        [Ts,XMGliq,test,testS,testSdry,CL,mgratio,NaKratio,WNACA,XLFG,XLFS,XLFP,WT,WK,WNA,S,Tsdry,PCA,QCA,XCP,XPP,XCS,XPS,XCG,XPG,Dn] = ...
            iterateMelt_QCA_Oct6(Frac,ii, Tb, Pb,XOW,test,testS,water,Co,CR,Do,P,mgratio,NaKratio,XMG,YMG,XMP,YMP,Ts,XMS,YMS,Tsdry,testSdry,Wo,WnS,Lnr,Wnr,DMP,DMS,DMG, Gar2SpTransition_pressure, Sp2PlagTransition_pressure,H2O_index);
        test_dFdP(ii) = iiii;
        iiii = iiii +1;
    else
        test_dFdP(ii) = 0;
    end
    
    testMODE(ii,:) = [XLFG QCA DCA PCA DG];
    testPCA(ii,1) = PCA;
    testPCA(ii,2) = XLFP(3).*DMP(3,1);
    Dn_Test(ii,:) = Dn;
    Do_out(ii,:) = Do;
    P_out(ii,:) = P;
    
    
    
    tempNaK_Plag =[ (XLFP(3)+XLFP(4))/(XLFP(3)+XLFP(4)+QCA)];
    testPboundary_Sp2Gar = (1.2  - .52*WNA+2.11*XMGliq)*10;
    testPboundary_Plag2Sp = (0.823  +1.138*WNA+0.165*XMGliq)*10;
    testPboundary_Plag2Sp_PlagSide = (0.823  +1.138*tempNaK_Plag+0.165*XMGliq)*10; %This is using the Mg\# from the spinel field
    
    if testPboundary_Plag2Sp_PlagSide < testPboundary_Plag2Sp
        testPboundary_Plag2Sp = testPboundary_Plag2Sp_PlagSide;
    end
    
    
    if Pb(ii) >= Gar2SpTransition_pressure
        if all(testPboundary_Sp2Gar >= min(Pb) & testPboundary_Sp2Gar <= max(Pb))==1
            tmp = abs(Pb-testPboundary_Sp2Gar);
            [idx Gar2SpTransition_index_temp] = min(tmp); %index of closest value
            Gar2SpTransition_pressure = Pb(Gar2SpTransition_index_temp) ;
            if (Pb(ii) - Gar2SpTransition_pressure) < 0
                Gar2SpTransition_pressure = Pb(ii);
            end
        else
            Gar2SpTransition_pressure=testPboundary_Sp2Gar;
        end
        
        
    else if Pb(ii) >= Sp2PlagTransition_pressure    %Spinel Field
            if all(testPboundary_Plag2Sp >= min(Pb) & testPboundary_Plag2Sp <= max(Pb))==1
                tmp = abs(Pb-testPboundary_Plag2Sp);
                [idx Sp2PlagTransition_index_temp] = min(tmp); %index of closest value
                Sp2PlagTransition_pressure = Pb(Sp2PlagTransition_index_temp);
                if (Pb(ii) - Sp2PlagTransition_pressure) < 0
                    Sp2PlagTransition_pressure = Pb(ii);
                end
            else
                Sp2PlagTransition_pressure = testPboundary_Plag2Sp;
            end
            
            
        end
        
    end
    
    if strcmp(forceSpPlag,'force Sp-Plag Boundary')==1
        Sp2PlagTransition_pressure = Sp2PlagTransition_forced;
    end
    
    S = S./sum(S);
    Sdry = Sdry./sum(Sdry);
    
    test(ii,:) = [Pb(ii) Ts(ii) Tb(ii) (1-XMGliq) WNA WT WK WNACA testPboundary_Plag2Sp testPboundary_Sp2Gar QCA Sp2PlagTransition_pressure Gar2SpTransition_pressure];
    testS(ii,:) = [S sum(S)];
    testSdry(ii,:) = Sdry;
    
    
    if isnan(Ts(ii)) == 1;
        if ii ~=1
            Ts(ii) = Ts(ii-1);
        else
            Ts(ii)=1500;
        end
    end
    
    clear CL
    
    if T(ii)<Ts(ii)
        % melting has not actually happened, but we still need to update the
        % masses related to melting
        
        % liquid composition is 0 in both major and trace elements
        BAI3(ii,:) = zeros(1,9);
        CL = 0.*CR;
        % residual solids do not change
        sXOW = sXOW;
        CS = CS;
        melton = 0;
        Frac = 0;
        rFrac = rFrac;
        if ii> 1
            ln(ii) =  Frac.*WnS(ii-1); %
            WnS(ii) = WnS(ii-1) - ln(ii); % the solid mass will stay the same (ln(ii) = 0)
            Ln(ii) = ln(ii) + Lnr(ii-1); % the melt mass will also stay the same (ln(ii) = 0)
            cumFrac =  (Wo-WnS(ii))./Wo;%1-WnS(ii);
        else
            ln(ii) = Frac.*Wo;
            WnS(ii) = Wo - ln(ii);
            Ln(ii) = ln(ii);
            cumFrac =  (Wo-Wo)./Wo;%1-WnS(ii);
        end
        CumFrac(ii) = cumFrac;
        rFac_out(ii) = rFrac;
        
        
    else
        % melting has now happened, update all the masses and calculate liquid
        % compositions
        
        
        IDFirstMelt=IDFirstMelt+1;
        if IDFirstMelt==1
            
            VdFdP = ones(1,length(Pb))*0.002;
            VdFdP(ii:end) = 0.002 + 0.023/Pb(ii)*(Pb(ii)-Pb(ii:end));
            
            if ismember(dFdPtype,'variable') ==1
                dFdP = VdFdP;
                
                % Reinitialize melt fraction
                MeltFrac = zeros(1,length(Pb));  %Array for melt fractions
                for i = 1:length(Pb)
                    if length(Pb) == 1
                        MeltFrac(i) = dFdP;
                    else if size(unique(Pb),2)==1
                            MeltFrac(i) = dFdP(i);
                        else
                            if i == 1
                                MeltFrac(i) = -(Pb(i+1)-Pb(i)).*dFdP(i);
                            else
                                if Pb(i)-Pb(i-1)==0
                                    fixPbx(i)=i;
                                    MeltFrac(i) = dFdP(i);
                                else
                                    MeltFrac(i) = -(Pb(i)-Pb(i-1)).*dFdP(i);
                                end
                            end
                        end
                    end
                end
                
                if exist('fixPbx')==1
                    MeltFrac(find(fixPbx,1)-1) = dFdP(find(fixPbx,1)-1);
                end
            end
        end
        
        
        
        melton = 1;
        Frac = MeltFrac(ii);
        
        if ii> 1
            
            ln(ii) =  Frac.*WnS(ii-1);
            WnS(ii) = WnS(ii-1) - ln(ii);
            Ln(ii) = ln(ii) + Lnr(ii-1);
            CL = CR.*Wnr(ii-1)./(Dn.*WnS(ii)+Ln(ii));
            
        else
            ln(ii) = Frac.*Wo;
            WnS(ii) = Wo - ln(ii);
            Ln(ii) = ln(ii);
            CL = Co.*Wo./(Dn.*WnS(ii)+Ln(ii));
            
        end
        rFrac  = Ln(ii)./(WnS(ii)+Ln(ii));
        
        % Calculate the melt composition
        TA = zeros(1,8);
        TA(1) = S(4); %quartz
        TA(3) = S(3); %plag
        TA(4) = S(1); %olivine
        TA(5) = S(2); % cpx
        
        % Renormalize Ol, Cpx, Plag, Qtz to 1.0
        %Qtz
        %Plag
        %Olivine
        %Cpx
        % Recalculate to mole components
        TA(1)=TA(1);     %Qtz
        TA(2)=TA(2)/4.;  %Albite
        TA(3)=TA(3)/4.;  %Anorthite
        TA(4)=TA(4)/2.;  %Olivine
        TA(5)=TA(5)/3.;  %Cpx
        TA(6)=TA(6)/1.5; %TiO2-bearing mineral
        TA(7)=TA(7)/4.;  %K-bearing mineral
        TA(8)=TA(8)/6.;  %
        TA = TA./sum(TA);
        
        %Recalculate albite & anorthite contents in plag component
        TPL=TA(2)+TA(3);  %plag
        XBRT=2.*WNACA/61.979;    %ratio1 of KG92
        XNRT=(1.-WNACA)/56.0794; %ratio2 of KG92
        XNA=XBRT/(XBRT+XNRT);
        TA(2)=XNA*(TPL+TA(5)); %albite
        TA(3)=TPL-TA(2); %anorthite
        
        
        GP(9)=TA(7);     %K2O (cation frac)
        GP(10)=TA(2);    %Na2O (cation frac)
        % Transform from molar to cation fraction
        GP(3)=2.*TA(3)+TA(2)+TA(7);  %Al2O3
        GP(8)=TA(5)+TA(3);   %CaO
        GP(2)=TA(6);  %TiO2
        GP(5)=(1-XMGliq)*(2.*TA(4)+TA(6)+TA(5)); %FeO
        GP(6)=XMGliq*(2.*TA(4)+TA(6)+TA(5));  %MgO
        GP(1)=TA(1)+TA(4)+2.*TA(5)+2.*TA(3)+3.*TA(7)+3.*TA(2); %SiO2
        GP(4)=0.0; %Cr2O3
        GP(7)=0.0; %MnO
        GP(11)=0.0;  %Fe2O3
        GP(12)=0.0;  %P2O5
        GP = GP./sum(GP);
        
        
        % Transform from cation fraction to mole %
        GN(1)=GP(1);        %SiO2
        GN(2)=GP(2);        %TiO2
        GN(3)=0.5*GP(3);    %Al2O3
        GN(4)=GP(4);        %Cr2O3
        GN(5)=GP(5);        %FeO
        GN(6)=GP(6);        %MgO
        GN(7)=GP(7);        %MnO
        GN(8)=GP(8);        %CaO
        GN(9)=0.5*GP(9);    %K2O
        GN(10)=0.5*GP(10);  %Na2O
        GN(11)=0.5*GP(11);  %Fe2O3
        GN(12)=0.5*GP(12);  %P2O5
        GN = 100.*GN./sum(GN);
        
        
        % Calculate composition in wt%
        CN = 100.*GN.*FK./sum(GN.*FK);
        
        % Add Ti, Cr, & K
        CN(2) = WT;   
        %Ti
        if Pb(ii)>=23
            CN(4) = XLFG(2);    %Cr
            %CN(10) = XLFG(4);
        else if Pb(ii) >= 9
                CN(4) = XLFS(2);    %Cr
            else
                CN(4) = XLFP(2);    %should be right now!
            end
        end
        
       
        CN(9) = WK;             %K
        CN = 100.*CN./sum(CN);
        
        
        
        % Transform to mole fraction
        GN = 100.*(CN./FK)./sum(CN./FK);
        % Recaculate molar oxides with addition of ilmenite and orthoclase
        % components
        % Things get realitivley more silica rich, more aluminum rich, and more
        % iron and magnesium rich.
        
        GN(1)=GN(1)+6.*GN(9);  %SiO2 = SiO2 + 6 Or
        GN(3)=GN(3)+GN(9);     %Al2O = Al2O3 + Or
        GN(5)=GN(5)+(1-XMGliq)*GN(2);  %FeO + (1-Mg#)*Ilm
        
        GN(6)=GN(6)+XMGliq*GN(2); %MgO + Mg# * Ilm
        %debugFe(ii,:) = [GN(5), 1-XMGliq,GN(2),GN(6) ];
        GN = 100.*GN./sum(GN);
        % debugFe(ii,:) = [GN];
        %Recalculate composition in wt%
        CN = GN.*FK;
        CN = 100.*CN./sum(CN);
        
        % Rearrange to remove MnO, Fe2O3, P2O5 from saved melt composition
        BIW = CN;
        BIW(7) = CN(8); %cao
        BIW(8) = CN(9); %k2o
        BIW(9) = CN(10); %na2o
        
        
        
        %%
        % Save instantaneous major and trace element compositions at each depth
        %[SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O]
        BAI3(ii,:)=100*BIW(1:9)./sum(BIW(1:9));
        CL_BAI3(ii,:) = CL;
        % Calculate cumulative melt fraction
        cumFrac = (Wo-WnS(ii))./Wo;%1-WnS(ii);
        CumFrac(ii) = cumFrac;
        IncFrac(ii) = Frac;
        rFrac_out(ii) = rFrac;
        

        
 
        % Calculate major and trace element solid residue compositions
        if ii> 1
            sXOW = (XOW.*Wnr(ii-1) - BAI3(ii,:).*Ln(ii))./WnS(ii);
            CS = (CR.*Wnr(ii-1) - CL.*Ln(ii))./WnS(ii);
            
        else
            sXOW = (XOW.*Wo - BAI3(ii,:).*Ln(ii))./WnS(ii);
            CS = (CR.*Wo - CL.*Ln(ii))./WnS(ii);
        end
        
        
        % % %         if any(sXOW<0)
        % % %             'solid residue ran out of an element, something wrong'
        % % %             return
        % % % %             find(sXOW<0)
        % % % %             sXOW(sXOW<0) = 0;
        % % % %             sXOW=sXOW./sum(sXOW)*100;
        % % %         end
        % % %
        
        
    end
    
    %%
    % No matter if melting happneed or not, calculate masses related to melt extraction
    if Ln(ii)/(WnS(ii)+Ln(ii)) <= rFrac_crit
        % If the mass of liquid is less than the critical melt fraction, don't extract anything:
        % Calculate bulk composition of mantle
        % need to check if things go negative SMB XX
        
        rFrac = rFrac;
        debugFe(ii,:) = [XOW(5:6) sXOW(5:6) BIW(5:6)];
        
        %don't extract any melt
        Qn(ii) = 1;
        Lnr(ii) = Qn(ii).*Ln(ii);
        Lnexp(ii) = (1-Qn(ii)).*Ln(ii);
        Wnr(ii) = WnS(ii)+Lnr(ii);
        
        % update residue (solid + liquid) composition, although it will not
        % change becuase there is no melt extraction
        if ii > 1
            XOW = ((XOW.*Wnr(ii-1))-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
        else
            XOW = ((XOW.*Wo)-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
        end
        CR = CR;
        %save extracted melt compositions
        BAI4(ii,:)=zeros(1,9);     % No extracted melt increment (set to 0)
        CL_OUT(ii,:)=0*CL;
        
    else % extract melt!
        
        if T(ii)>=Ts(ii)
            if ii> 1
                Qn(ii) = (rFrac_crit.*WnS(ii))./(ln(ii).*(1-rFrac_crit)+Lnr(ii-1).*(1-rFrac_crit));
                Lnr(ii) = Qn(ii).*Ln(ii);
                Lnexp(ii) = (1-Qn(ii)).*Ln(ii);
                Wnr(ii) = WnS(ii)+Lnr(ii);
                % update residue (solid + liquid) major and trace element compositions
                XOW = ((XOW.*Wnr(ii-1))-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
                CR =  ((CR.*Wnr(ii-1))-(CL.*Lnexp(ii)))./Wnr(ii);
                
            else
                Qn(ii) = (rFrac_crit.*WnS(ii))./(ln(ii).*(1-rFrac_crit)+Lnr(1).*(1-rFrac_crit));
                Lnr(ii) = Qn(ii).*Ln(ii);
                Lnexp(ii) = (1-Qn(ii)).*Ln(ii);
                Wnr(ii) = WnS(ii)+Lnr(ii);
                % update residue (solid + liquid) major and trace element compositions
                XOW = ((XOW.*Wo)-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
                CR =  ((CR.*Wo)-(CL.*Lnexp(ii)))./Wnr(ii);
            end
            
            %compare to 2D melt matrix, assumes Lnexp is the same and melt
            %composition are the same if the temperature difference is
            %comparable to the temeprature above the solidus.  This is not
            %calcualted correctly and is only a back of the envelope check.
            Tb_IT_compare = Tb(ii) - IT_extended(ii,:);
            TaboveTs= T(ii) - Ts(ii);
            
            for gg = 1:size(IT_extended,2)
                if TaboveTs >=Tb_IT_compare(gg)
                    Lnexp_IT(ii,gg) = Lnexp(ii);
                else
                    Lnexp_IT(ii,gg) = 0;
                end
            end
            
            
            %save extracted melt compositions
            BAI4(ii,:)=100*BIW(1:9)./sum(BIW(1:9));  %  Extracted melt increment
            CL_OUT(ii,:)=  CL;
            rFrac = rFrac_crit;
            debugFe(ii,:) = [XOW(5:6) sXOW(5:6) BIW(5:6)];
            
        else
            
            rFrac = rFrac;
            debugFe(ii,:) = [XOW(5:6) sXOW(5:6) BIW(5:6)];
            
            %don't extract any melt
            Qn(ii) = 1;
            Lnr(ii) = Qn(ii).*Ln(ii);
            Lnexp(ii) = (1-Qn(ii)).*Ln(ii);
            Wnr(ii) = WnS(ii)+Lnr(ii);
            
            % update residue (solid + liquid) composition, although it will not
            % change becuase there is no melt extraction
            if ii > 1
                XOW = ((XOW.*Wnr(ii-1))-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
            else
                XOW = ((XOW.*Wo)-(BAI3(ii,:).*Lnexp(ii)))./Wnr(ii);
            end
            CR = CR;
            %save extracted melt compositions
            BAI4(ii,:)=zeros(1,9);     % No extracted melt increment (set to 0)
            CL_OUT(ii,:)=0*CL;
        end
        
    end
    
    %% update mantle modes based on melting reaction
    
    for i=1:6
        XMP(i)=(XMP(i)-Frac*YMP(i));
        XMS(i)=(XMS(i)-Frac*YMS(i));
        XMG(i) =(XMG(i)-Frac*YMG(i));
        
        if XMP(i) <= 0; XMP(i) = 0; end
        if XMS(i) <= 0; XMS(i) = 0;  end
        if XMG(i) <= 0; XMG(i) = 0; end
    end
    
    XMP = XMP./sum(XMP);
    XMS = XMS./sum(XMS);
    XMG = XMG./sum(XMG);
    
    
    %% Convert remaining Gar to Sp at P = 9 kbar
    
    if Pb(ii) == Gar2SpTransition_pressure
        Gar2SpTransition_index = ii;
        
        
        
        %for CaMg2Al2Si3O12
        %         XMS(6) = 0;
        %         XMS(4) = 0;
        %         XMS(5) = (3/8)*XMG(6);
        %         XMS(1) = XMG(1) + (1/2)*XMG(6);
        %         XMS(2) = XMG(2) + (9/8)*XMG(6);
        %         XMS(3) = XMG(3) - XMG(6);
        
        %for (Ca1/7Mg6/7)3Al2Si3O12
        XMS(6) = 0;
        XMS(4) = 0;
        XMS(5) = (3/8)*XMG(6);
        XMS(1) = XMG(1) + (3/14)*XMG(6);
        XMS(2) = XMG(2) + (79/56)*XMG(6);
        XMS(3) = XMG(3) - XMG(6);
        
        % not an ideal solution, but neccessary if PCA equation is not smooth
        %Fe-Mg Kds
        XOlK = 0.31509+0.0062362*(Pb(ii)/10);
        XCpxK = 0.27652+0.021516*(Pb(ii)/10);
        XOpxK = 0.27262+0.042031*(Pb(ii)/10);
        %XSpK = 1.05210-0.311430*(Pb(ii)/10);
        XSpK = 0.6;    %Spinel Fe-Mg Kd
        if ii>1 && ii~=max(ii)
            XMGliq = BG15_Mgnum_MassBalance(Frac, Wnr(ii-1),WnS(ii), XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
            [PCA] = calculatePCA(Pb(ii+1), XMGliq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
        else
            XMGliq = BG15_Mgnum_MassBalance(Frac, Wo,Wo - (Frac.*Wo), XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
            [PCA] = calculatePCA(Pb(ii), XMGliq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
        end
        %PCA = 24.49-0.469*Pb(ii)-12.92*(1-XMGliq); % needed below to recalc mode
        
        
    end
    
    % Convert remaining Sp to Plag at P = 9 kbar
    if Pb(ii) == Sp2PlagTransition_pressure
        Sp2PlagTransition_index = ii;
        XMP(4) = 1.5*XMS(5);
        XMP(5) = 0;
        XMP(1) = XMS(1) - XMS(5);
        XMP(2) = XMS(2) - XMS(5);
        XMP(3) = XMS(3) + 1.5*XMS(5);
        %'converted remaining spinel to plag'
    end
    
    
    % update opx - cpx proportions according to the lever rule
    %if 3 >  30   % uncomment to prevent this loop from occuring
    clear Pyx
    Pyx = XMS(1) + XMS(2);
    if Pyx>0
        Cabulkfactor = 1./Pyx;
        XMS(1) = Pyx*(Cabulkfactor*XOW(7) - CaOOPX)./(PCA - CaOOPX);
        XMS(2) = Pyx-XMS(1);
        if XMS(1) <= 0; XMS(1) = 0;  XMS(2) = Pyx; end
        if XMS(2) <= 0; XMS(2) = 0;  XMS(1) = Pyx;end
        XMS = XMS./sum(XMS);
        % Calculate initial pyroxene modes in garnet field based on P-dependent
        %   cpx-opx join
        testXMG(ii,:) = XMG(1:2) ;
        clear Pyx PPX
        Pyx = XMG(1) + XMG(2);
        Cabulkfactor = 1./Pyx;
        PPX = Cabulkfactor*(XOW(7) - 4.96*XMG(6));
        [PCA] = calculatePCA(Pb(ii), XMGliq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
        % if Pb(ii)>=40
        %     Pbhold = 40;
        %     PCA = 20.167-11.6383*(1-MgLiqEstimate)-0.23985*(Pbhold);
        % else
        %     PCA = 20.167-11.6383*(1-XMGliq) - 0.23985*Pb(ii);
        % end
        % %PCA = 20.167-11.6383*(1-XMGliq) - 0.23985*Pb(ii);
        
        XMG(1) = Pyx*(PPX - CaOOPX)./(PCA - CaOOPX);
        XMG(2) = Pyx-XMG(1);
        
        if XMG(1) <= 0; XMG(1) = 0; XMG(2) = Pyx;  end
        if XMG(2) <= 0; XMG(2) = 0; XMG(1) = Pyx; end
        
        XMG = XMG./sum(XMG);
        testXMG(ii,:) = XMG(1:2) - testXMG(ii,:);
        Cabulk_mode(ii+1) = PPX;
        
    end
    
    % If plag exhausted correct melting reaction
    if XMP(4) == 0 && Pb(ii) <= 9
        YMP(1) = 0.66;
        YMP(2) = -0.02;
        YMP(3) = 0.36;
        YMP(4) = 0;
        YMP(5) = 0;
        YMP(6) = 0;
    end
    
    % Save mantle mode
    if Pb(ii) > Gar2SpTransition_pressure
        XM_OUT(ii,:) = XMG;
        testMODEtype(ii) = 1;
    else if Pb(ii) > Sp2PlagTransition_pressure
            XM_OUT(ii,:) = XMS;
            testMODEtype(ii) = 2;
        else
            XM_OUT(ii,:) = XMP;
            testMODEtype(ii) = 3;
        end
    end
    
    
    
    
    CheckPhaseExistance = nnz(~XM_OUT(ii,:));
    if CheckPhaseExistance==5
        'Wow, all phases but one were exhuasted'
        [aaa,bbbb]=find(XM_OUT(ii,:)>0);
        YMS(:)=0;
        YMP(:)=0;
        YMG(:)=0;
        YMS(bbbb)=1;
        YMP(bbbb)=1;
        YMG(bbbb)=1;
    end
    
    
    
    
    XOW_OUT(ii,:) = XOW;
    sXOW_OUT(ii,:) = sXOW;
    CS_OUT(ii,:) = CS;
    CR_OUT(ii,:) = CR;
    % distribute trace elements among the remaining solid phases
    % calculate trace element melting
    if Pb(ii) >= Gar2SpTransition_pressure
        Do = Dmatrix*XMG';
    else if Pb(ii) >= Sp2PlagTransition_pressure
            Do = Dmatrix*XMS';
        else
            Do = Dmatrix*XMP';
        end
    end
    
    
    CS_minerals = bsxfun(@times, (bsxfun(@rdivide, Dmatrix,Do)),CS);
    CS_minerals_OUT(ii,:,:) = CS_minerals;
    %for calculating massbalance on trace elements
    %(CR_OUT_1(25,:).*Wnr(25)) - (CL_BAI3_1(25,:).*Ln_1(25) + CS_OUT_1(25,:).*WnS(25))
    
    
    
    
    
end


if ~exist('Gar2SpTransition_index')
    Gar2SpTransition_index=ii;
end
if ~exist('Sp2PlagTransition_index')
    Sp2PlagTransition_index=ii;
end
Fc = max(CumFrac);
iFc = find(CumFrac == Fc);




% if length(iFc) > 1
%     CumFrac(iFc(2:end)) = 0;
% end

% size(CL_BAI3)
% size(CR_OUT)
% size(CS_minerals_OUT)
%
% size(helpDEBUG)
% size(XM_OUT)
% size(XOW_OUT)
% size(BAI3)
% size(BAI4)
% size(Ts)
% size(T)
% size(CumFrac)
% size(IncFrac)
% size(Fc)
% size(sXOW_OUT)
% size(mgratio)
% size(NaKratio)
% size(test)
% size(testS)
% size(debugFe)
% size(CL_OUT)
% size(CS_OUT)
% size(CR_OUT)
% size(CL_BAI3)
% size(CS_minerals_OUT)
% size(testSdry)
% size(test_dFdP)
% size(testXMG)
% size(testMODE)
% size(testPCA)
% size(Cabulk_mode)
% size(MeltFrac)
% size(Do_out)
% size(P_out)
% size(Dn_Test)
% ln
% ,ln
% WnS
% Ln, Qn, Lnr, Lnexp,
% Wnr,testMODEtype,initialBulkComp,KDs,SolidResidueDensity, LiquidDensity,Gar2SpTransition_index, Sp2PlagTransition_index
CS_minerals_OUT = CS_minerals_OUT(:,:,1);

end


%% Additional Subroutines

function [RT4BIS] = BG15_Mgnum_MassBalance(Frac, Wnminus1R,Wns, XM,YM,XOlK,XCpxK,XOpxK,XSpK,XGarK,XO)
%==========================================================================
% This is a modified version of the rtbis routine found in
% Numerical Recipes.  Its purpose is to be an elegant root finder

% This version is modified to estimate the mg# of melt from mantle
% assumed to consist of pl or sp lherz as function of F (extent of melt)
% and depends on ol/melt KD and cpx/melt, opx/melt, and sp/melt KDs
% and phase proportions in mantle and on bulk mgo and feo
%
% FL is melt fraction (WT), X1 and X2 are guesses of mgnum of melt, XACC is
% tolerance, XOlK, XCpxK, XOpxK, XSpK are the ol/melt cpx/melt, opx/melt,
% and sp/melt KDs, repsectively, XM(5) is vector of phase proportions (WT)
% MMGNUM is function that solves for melt mgnum
%
% This version of the RT4BIS subroutine from Kinzler & Grove 92 takes XO,
%   XM, YM, and FL in weight fractions.  Note that the order of elements in
%   XO is [K,Na,Ca,Fe,Mg,Ti,Si,Al] as in the original KG92 routines.
% An important component are the assumptions that go into the ol,px,&sp
%   composition in FUNC1-FUNC4 in the MMGNUM subroutine.
% This routine is valid in both the plag & spinel fields because plag will
%   not affect the Mg-Fe partitioning (i.e., XMM(4) does not come into the
%   final calculation of Mg#)
% M. Behn & T. Grove 11/13/12
% M. Behn adjusted Ca/Mg+Fe for cpx from 25 (KG92) to 40 (K97) 11/14
%==========================================================================
XACC = 1e-5;
X1 = 0.85;
X2 = 0.50;
%X2 = 0.68;
% X1 = 0.80;
% X2 = 0.50;
%     G1 = 0.80;
%     G2 = 0.50;

JMAX=40;
X21diff = 0.001;
X21diff_fine = 100;
%JMAX=2000;

% Use melt fraction (FL), initial phase proportions (XM) and melt
% reaction (YM) to estimate phase proportions at end of melting (XMN) (all
% in wt. units)

for i=1:6
    XMN(i)=(XM(i)-Frac*YM(i));
    if XMN(i) <= 0;
        XMN(i) = 0;
    end
end


XMN = XMN./sum(XMN);
XFEOM = XO(5);
XMGOM = XO(6);
F=MMGNUM_MassBalance(X1,Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK);



if F < 0
    RT4BIS=X1;
    DX=X2-X1;
else
    RT4BIS=X2;
    DX=X1-X2;
end


for J = 1:JMAX
    DX=DX*.5;
    XMID=RT4BIS+DX;
    FMID=MMGNUM_MassBalance(XMID,Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK);
    
    if FMID < 0;
        RT4BIS=XMID;
    end
    
    if abs(FMID) < XACC ||  FMID == 0;
        RT4BIS=XMID; %SB: put this here to prevent errors occuring when the F is initially < 0.
        return;
    end
end

if J == JMAX
    
    
    
    
    triedx = X2:X21diff:X1;
    RT4BIS = zeros(1,length(triedx));
    for ii=1:length(triedx)
        RT4BIS(ii) = MMGNUM_MassBalance(triedx(ii),Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK);
    end
    %val = 0; %value to find
    %tmp = abs(RT4BIS-triedx-val)
    %[idx idy] = min(abs(RT4BIS-triedx)); %index of closest value
    [idx idy] = min(abs(RT4BIS)); %index of closest value
    %closest = RT4BIS(idy) %closest value
    %
    %                         figure()
    %                 hold on
    %                 triedx
    %                 RT4BIS
    %                 plot(triedx, RT4BIS,'*')
    %                 hline(0)
    %                 hline(idx)
    %                 xlabel('triedx')
    %                 ylabel('RTBIS, where RTBIS=RTBIS-MGNUM');
    
    
    if idy == 1 || idy == length(triedx)
        
        
        % 'TESTTEST'
        %
        %
        % Wns
        % Wnminus1R
        % XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK
        % idy
        % length(triedx)
        % YM
        % return
        triedx = X2:X21diff/X21diff_fine:X1;
        RT4BIS = zeros(1,length(triedx));
        for ii=1:length(triedx)
            RT4BIS(ii) = MMGNUM_MassBalance(triedx(ii),Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK);
        end
        [idx idy] = min(abs(RT4BIS)); %index of closest value
        RT4BIS = triedx(idy);
        
    else
        triedx = triedx(idy-1):X21diff/X21diff_fine:triedx(idy+1);
        RT4BIS = zeros(1,length(triedx));
        for ii=1:length(triedx)
            RT4BIS(ii) = MMGNUM_MassBalance(triedx(ii),Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK);
        end
        %val = 0; %value to find
        %tmp = abs(RT4BIS-triedx-val)
        [idx idy] = min(abs(RT4BIS)); %index of closest value
        %closest = RT4BIS(idy) %closest value
        RT4BIS = triedx(idy);
    end
end
end

function [FUNC] = MMGNUM_MassBalance(MGNUM,Wnminus1R,Wns, XMGOM,XFEOM,XMN,XOlK,XCpxK,XOpxK,XSpK,XGarK)

Mg_nminus1s = XMGOM;
Fe_nminus1s = XFEOM;
mantlemode = XMN;
Mgnum_nl_moles = MGNUM;


% convert Mg#melt guess to wt% for comparison at the end
p1 = 0.1091;
p2 = -0.089921;
p3 = 0.19553;
p4 = 0.22155;
p5 = 0.56373;
p6 = -5.3689e-05;
x = Mgnum_nl_moles;
Mgnum_nl_wtpercent = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6;

KDs = [XCpxK XOpxK XOlK  0 XSpK XGarK];

%define the composition of the minerals using stoichiometry (setting FeO = 0 for now) to eventually
%calcualte the wt% of MgO and wt% of FeO in the solid

%[1-SiO2,2-TiO2,3-Al2O3,4-Cr2O3,5-FeO,6-MgO,7-CaO,8-K2O,9-Na2O]
atomicoxideweights = [60.084 79.866 101.961 151.9902 71.844 40.304 56.077 94.196 61.979];

stoichoimetricMinerals = [...
    2.0	1.0	1.0	2.0	0.0	3.0 %% SiO2
    0.0	0.0	0.0	0.0	0.0	0.0 %% TiO2
    0.0	0.0	0.0	1.0	1.0	1.0 %% Al2O3 %fixed 10/3/17
    0.0	0.0	0.0	0.0	0.0	0.0 %% Cr2O3
    0.0	0.0	0.0	0.0	0.0	0.0 %% FeO : for right now assume no FeO
    1.6	1.0	2.0	0.0	1.0	3.0 %% MgO
    0.4	0.0	0.0	1.0	0.0	0.0 %% CaO %fixed plag 8/18/18
    0.0	0.0	0.0	0.0	0.0	0.0 %% K2O
    0.0	0.0	0.0	0.0	0.0	0.0]; %% Na2O

% convert mineral formulas to atomic percent
sumatomicpercentMinerals=sum(stoichoimetricMinerals);
atomicpercentMinerals = bsxfun(@rdivide, stoichoimetricMinerals,sumatomicpercentMinerals).*100;

% using the input Mg# melt guess in moles, calculate the Mg# of each mineral phase
Mgnum_ns = Mgnum_nl_moles./(Mgnum_nl_moles+KDs.*(1-Mgnum_nl_moles));

% update atomicpercentMinerals with the correct partitioning of FeO and MgO
atomicpercentMinerals(5,:) = atomicpercentMinerals(6,:).*(1-Mgnum_ns); % FeO
atomicpercentMinerals(6,:) = atomicpercentMinerals(6,:).*Mgnum_ns; % MgO

%convert atomic percents to wt%
atomicoxideweightstranspose = atomicoxideweights';
wtpercentminerals = bsxfun(@times, atomicpercentMinerals,atomicoxideweightstranspose);
sumwtpercentminerals = sum(wtpercentminerals);
wtpercentminerals = bsxfun(@rdivide, wtpercentminerals,sumwtpercentminerals).*100;

%sum up the MgO and FeO in each of the mineral phases to calcualte the bulk MgO and FeO in the solid in wt%
Mg_ns = mantlemode*wtpercentminerals(6,:)';
Fe_ns = mantlemode*wtpercentminerals(5,:)';

% % %backtrack by returing to moles to check to make sure I get the same answer
% % %as the previous script
% % bulk_ns_wtpercent = mantlemode*wtpercentminerals';
% % bulk_ns_moles = bulk_ns_wtpercent./atomicoxideweights;
% % bulk_ns_moles = bulk_ns_moles./sum(bulk_ns_moles).*100;
% % Mg_ns_moles = bulk_ns_moles(6);
% % Fe_ns_moles = bulk_ns_moles(5);
% % Mgnum_nl_RS_moles = (Mg_nminus1s - (1 - FLM).*Mg_ns_moles)./((Mg_nminus1s+Fe_nminus1s) - (1-FLM).*((Mg_ns_moles+Fe_ns_moles)))

%calculate the Mg# of the melt
Mgnum_nl_RS_wtpercent = ((Mg_nminus1s*Wnminus1R) - (Mg_ns*Wns))./(((Mg_nminus1s+Fe_nminus1s)*Wnminus1R) - (Wns*(Mg_ns+Fe_ns)));

%calculate the residual between the guess and the RHS of the equation, in wt%
FUNC=Mgnum_nl_RS_wtpercent-Mgnum_nl_wtpercent;
end

%%

function [XMG,XMS,XMP] = calculateInitialModes(ZO, P,KD,Gar2SpTransition_pressure,Sp2PlagTransition_pressure)

P_GST = Gar2SpTransition_pressure;
P_SPT=  Sp2PlagTransition_pressure;

% 
% %script to calculate a mode given a matrix of bulk compositions
% %Mantle Composition [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;
% % 1 - HZ prim (KG93)
% %ZO(1,:) =  [46.16,0.182,4.08,0.47,7.57,37.95,3.22,0.03,0.333];  %prim (KG93)% %% grove
% ZO(1,:) =  [46.161,0.182,4.08,0.47,7.573,37.946,3.224,0.032,0.333]; %% BG
% % 2 - HZ dep 1 (KG93)
% ZO(2,:) = [46.331,0.174,3.935,0.47,7.659,38.452,3.193,0.0074,0.281];
% % 3 - WH DMM (EPSL05)
% ZO(3,:) = [44.71,0.13,3.98,0.57,8.18,38.73,3.17,0.006,0.28];
% % 4 - HZ dep 2
% ZO(4,:) = [46.313,0.161,3.639,0.47,7.670,38.890,3.101,0.003,0.223];
% % 5 - WH DMM (EPSL05) - 1% melt (calculated from BG15)
% ZO(5,:) = [44.812,0.124,3.864,0.471,8.229,39.124,3.136,1.6849e-06,0.24];
% % 6 - WH DMM (EPSL05) - 2% melt (calculated from BG15)
% ZO(6,:) = [44.758,0.118,3.722,0.474,8.252,39.399,3.082,1.5208e-09,0.195];
% % 7 - MM3 (Baker & Stolper 1994; Baker et al., 1995)
% ZO(7,:) = [45.5,0.11,3.98,0.68,7.18,38.3,3.57,0,0.31];
% % % 8 - Kushiro 1996
% % ZO(8,:) = [43.7,0.25,2.75,.28,1.38*.9+8.81,37.22,3.26,.14,.33];
% % % 9 - HZ in CMAS
% % ZO(9,:) = [46.381,0.0,4.099,0,0,42.395,3.559,0,0];

% normalize the bulk composition
for i = 1:size(ZO)
 ZO(i,:)=100*ZO(i,:)./sum(ZO(i,:));
end

ZO(:,10)=ZO(:,5) + ZO(:,6); 

MgNumBulk = ((ZO(:,6)/40.3044))./(((ZO(:,6)./40.3044)+((ZO(:,5)./71.8464))));
MgLiqEstimate = 1 ./ (1 + (1/KD).*((1-MgNumBulk)./MgNumBulk));

%matrix containing the compositions of garnet stability field in order:
%olivine, opx, cpx, garnet

% compositionmatrixmaster = [40.3	0	0.21	0.14	9.48	49.59	0.3	0	0
% 54.74	0.07	5.48	0.64	5.16	31.47	2.42	0	0.23
% 54.27	0.07	5.56	0.61	5.08	26.47	7.37	0	0.59
% 42.48	0.33	22.29	1.66	5.46	22.98	4.96	0	0.02]; %comps from 40.07 Walter1998

compositionmatrixmaster = [...
40.3	0	0.13	0.08	10.32	49.04	0.21	0	0
55.61	0.14	3.94	0.3	6.14	32.14	1.93	0	0.27
55.4	0.2	4.09	0.39	5.33	24.93	9.59	0	0.83
42.47	0.52	22.34	1.32	6.45	21.99	4.71	0	0.04]; %comps from 40.02 Walter1998

% 
% compositionmatrixmaster = [40.3	0	0.21	0.14	9.48	49.59	0.3	0	0
% 56.00	0.07	4	0.64	5.16	31.47	2.42	0	0.23
% 56.00	0.07	4	0.61	5.08	26.47	7.37	0	0.59
% 42.48	0.33	22.29	1.66	5.46	22.98	4.96	0	0.02]; %made up comps

% 
% compositionmatrixmaster = [40.3	0	0.21	0.14	9.48	49.59	0.3	0	0
% 57.3700	0.07	2.63	0.64	5.16	31.47	2.42	0	0.23
% 57.3700	0.07	2.63	0.61	5.08	26.47	7.37	0	0.59
% 42.48	0.33	22.29	1.66	5.46	22.98	4.96	0	0.02]; %made up comps


% compositionmatrixmaster = [40.24	0	0.21	0.14	10.23	49.45	0.25	0	0
% 53.28	0.18	8.19	0.6	6.01	30.23	2.29	0	0.21
% 52.21	0.25	8.7	0.74	5.13	22.32	10.65	0	0.71
% 42.42	0.53	22.2	1.53	6.22	21.81	5.43	0	0.03]; %comps from 30.05 Walter1998





% normalize the mineral compositions? I think this makes sense to do
for i = 1:size(compositionmatrixmaster)
 compositionmatrixmaster(i,:)=100*compositionmatrixmaster(i,:)./sum(compositionmatrixmaster(i,:));
end
compositionmatrixmaster(:,10) = compositionmatrixmaster(:,5)+compositionmatrixmaster(:,6);


% reorder: cpx opx ol garnet
compositionmatrix(1,:) = compositionmatrixmaster(3,:);
compositionmatrix(2,:) = compositionmatrixmaster(2,:);
compositionmatrix(3,:) = compositionmatrixmaster(1,:);
compositionmatrix(4,:) = compositionmatrixmaster(4,:);


%select elements
% %CMAS
 %keep = [1 3 6 7];
 
%CMASF
%keep = [1 3 5 6 7];

%MASF
%keep = [1 3 5 6];


%C(M+F)AS
keep = [1 3 7 10];


compositionmatrix_CMAS = compositionmatrix(:,keep);

ZO_CMAS = ZO(:,keep);
%set normalize to 2 if you don't want to normalize, set it to 1 if you do
%want to normalize 
normalize = 1;

%renomalize the phase compositional matrix with respect to CMAS
if normalize == 1
    for i = 1:size(compositionmatrix_CMAS,1)
        compositionmatrix_CMAS(i,:) = 100*compositionmatrix_CMAS(i,:)./sum(compositionmatrix_CMAS(i,:));
    end
end

for i = 1:size(ZO)
%     alldata = regstats(ZO(i,:)',[cpx' opx' olivine' spinel'],'linear');
%     B = alldata.beta
%     sqrt(diag(alldata.covb))
%     alldata.rsquare
%     alldata.tstat.dfe
if normalize == 1 % normalize the bulk composition with respect to CMAS
    ZO_CMAS(i,:)=100*ZO_CMAS(i,:)./sum(ZO_CMAS(i,:));
end



% do some linear algebra i.e. mass balance and calculate an initial
% garnet mode
XMG(i,:) = compositionmatrix_CMAS'\ZO_CMAS(i,:)';
     
% if the mass balance results in any negative compostiions, redo by
% removing that mineral 
 if any(XMG(i,:)<0) == 1
     rowsremove = find(XMG(i,:)<=0);
     rowskeep = find(XMG(i,:)>0);
     compositionmatrix_CMAS_noopx =  compositionmatrix_CMAS;
     compositionmatrix_CMAS_noopx(rowsremove,:) = [];
     XMG_hold = compositionmatrix_CMAS_noopx'\ZO_CMAS(i,:)';
     XMG(i,rowskeep) =  XMG_hold;
     XMG(i,rowsremove) =  0;
 end
     
 %renormalize
 XMG(i,:) = XMG(i,:)./sum(XMG(i,:));
    
end


%enter the garnet mode in the 6 column matrix used in the main code
XMG = [XMG(:,1:3) zeros(size(XMG,1),2) XMG(:,4)];


%check if 1 or 2 pyroxenes
CaOOPX = 2;
PCA_GS = calculatePCA(P_GST, MgLiqEstimate, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
PCA_SP = calculatePCA(P_SPT, MgLiqEstimate, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);

%%
%if final P is in the garnet field, calculate the equlibrium pyroxene mode
if P>=P_GST
    %calcualte opx-cpx proportion given P
    PCA= calculatePCA(P, MgLiqEstimate, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
    
    clear Pyx
    Pyx = XMG(:,1) + XMG(:,2);
    Cabulkfactor = 1./Pyx;
    PPX = Cabulkfactor.*(ZO(:,7) - 4.96.*XMG(:,6));
    XMG(:,1) = Pyx.*((PPX - CaOOPX)./(PCA - CaOOPX));
    XMG(:,2) = Pyx-XMG(:,1);
    for i = 1:size(XMG,1)    
        if XMG(i,1) <= 0.; XMG(i,1) = 0.0; XMG(i,2) = Pyx(i);  end
        if XMG(i,2) <= 0.; XMG(i,2) = 0.0; XMG(i,1) = Pyx(i); end
    end
    
end


%%
% calculate spinel mode at transition

%first calculate the equlibrium pyroxene mode at the transition
clear Pyx PPX Cabulkfactor
Pyx = XMG(:,1) + XMG(:,2);
Cabulkfactor = 1./Pyx;
PPX = Cabulkfactor.*(ZO(:,7) - 4.96.*XMG(:,6));
XMG_GST(:,1) = Pyx.*((PPX - CaOOPX)./(PCA_GS - CaOOPX));
XMG_GST(:,2) = Pyx-XMG_GST(:,1);  
for i = 1:size(XMG,1)    
    if XMG_GST(i,1) <= 0.; XMG_GST(i,1) = 0.0; XMG_GST(i,2) = Pyx(i);  end
    if XMG_GST(i,2) <= 0.; XMG_GST(i,2) = 0.0; XMG_GST(i,1) = Pyx(i); end
end

%now do the spinel to garnet phase transition calculation
for i = 1:size(XMG,1)
    XMS(i,6) = 0;
    XMS(i,4) = 0;
    XMS(i,5) = (3/8)*XMG(i,6);
    XMS(i,1) = XMG_GST(i,1) + (3/14)*XMG(i,6);
    XMS(i,2) = XMG_GST(i,2) + (79/56)*XMG(i,6);
    XMS(i,3) = XMG(i,3) - XMG(i,6);
    XMS(i,:) = XMS(i,:)./sum(XMS(i,:));
end
        
%%
%if final P is in the spinel field, calculate the equlibrium pyroxene mode
if  P<P_GST && P >= P_SPT
    PCA= calculatePCA(P, MgLiqEstimate, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);        
    clear Pyx PPX Cabulkfactor
    Pyx = XMS(:,1) + XMS(:,2);
    Cabulkfactor = 1./Pyx;
    PPX = Cabulkfactor.*(ZO(:,7));
    XMS(:,1) = Pyx.*((PPX - CaOOPX)./(PCA - CaOOPX));
    XMS(:,2) = Pyx-XMS(:,1); 
     for i = 1:size(XMS,1)    
        if XMS(i,1) <= 0.; XMS(i,1) = 0.0; XMS(i,2) = Pyx(i);  end
        if XMS(i,2) <= 0.; XMS(i,2) = 0.0; XMS(i,1) = Pyx(i); end
    end   
end


%%
% calculate plag mode at transition

%first calculate equlibrum pyroxene mode at the transition 
clear Pyx PPX Cabulkfactor
Pyx = XMS(:,1) + XMS(:,2);
Cabulkfactor = 1./Pyx;
PPX = Cabulkfactor.*(ZO(:,7));
XMS_SP(:,1) = Pyx.*((PPX - CaOOPX)./(PCA_SP - CaOOPX));
XMS_SP(:,2) = Pyx-XMS_SP(:,1); 
 for i = 1:size(XMS,1)    
    if XMS_SP(i,1) <= 0.; XMS_SP(i,1) = 0.0; XMS_SP(i,2) = Pyx(i);  end
    if XMS_SP(i,2) <= 0.; XMS_SP(i,2) = 0.0; XMS_SP(i,1) = Pyx(i); end
 end   

%now calculate the transition reaction 
for i = 1:size(XMS,1)
XMP(i,4) = 1.5*XMS(i,5);
XMP(i,5) = 0;
XMP(i,1) = XMS_SP(i,1) - XMS(i,5);
XMP(i,2) = XMS_SP(i,2) - XMS(i,5);
XMP(i,3) = XMS(i,3) + 1.5*XMS(i,5);
XMP(i,6) = 0;
XMP(i,:) = XMP(i,:)./sum(XMP(i,:));
end
      

end

%%

function [PCA, DCA] = calculatePCA( Pb, MgLiq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure)
Pbhold = 30;


%PCA is CaO content in phases
%DCA is CaO D

if Pb>=Pbhold
    %     Pbhold = 23; %% SB changed from 27
    %     PCA = 24.49-0.469*Pbhold-12.92*(1-MgLiqEstimate);   
    PCA = 20.167-11.6383*(1-MgLiq)-0.23985*(Pbhold);
    DCA = 1.688 - 0.3486*(1-MgLiq) -0.0189*Pbhold;
    
else if Pb>=Gar2SpTransition_pressure
    PCA = 20.167-11.6383*(1-MgLiq)-0.23985*(Pb);
    DCA = 1.688 - 0.3486*(1-MgLiq) - 0.0189*Pb;
else if Pb>= Sp2PlagTransition_pressure % spinel melting
    PCA = 24.49-0.469*Pb-12.92*(1-MgLiq);
    DCA = 1.6 - 0.0196*Pb+0.071*(1-MgLiq);
else  % plag melting -- it's the same for plag and spinel 
    PCA = 24.49-0.469*Pb-12.92*(1-MgLiq);
    DCA = 1.6 - 0.0196*Pb+0.071*(1-MgLiq);
end
end
end

end

%%

function [Ts,XMGliq,test,testS,testSdry,CL,mgratio,NaKratio,WNACA,XLFG,XLFS,XLFP,WT,WK,WNA,S,Tsdry,PCA,QCA,XCP,XPP,XCS,XPS,XCG,XPG,Dn] = ...
iterateMelt_QCA_Oct6(Frac,ii, Tb, Pb,XOW,test,testS,water,Co,CR,Do,P,mgratio,NaKratio,XMG,YMG,XMP,YMP,Ts,XMS,YMS,Tsdry,testSdry,Wo,WnS,Lnr,Wnr,DMP,DMS,DMG,Gar2SpTransition_pressure, Sp2PlagTransition_pressure,H2O_index)

    if ii> 1
    Wnriiminus1   = Wnr(ii-1);
    ln_temp =  Frac.*WnS(ii-1);
    WnS_temp = WnS(ii-1) - ln_temp;        
    Ln_temp = ln_temp + Lnr(ii-1); 
    %Dn = (Do - P.*(Ln_temp))./(1-(Ln_temp));
    Dn = (Do - P.*(Frac))./(1-(Frac));
    CL = CR.*Wnr(ii-1)./(Dn.*WnS_temp+Ln_temp);        
    else
    Wnriiminus1   = Wo;
    ln_temp = Frac.*Wo;
    WnS_temp = Wo - ln_temp;       
    Ln_temp = ln_temp; 
    %Dn = (Do - P.*(Ln_temp))./(1-(Ln_temp));
    Dn = (Do - P.*(Frac))./(1-(Frac));
    CL = Co.*Wo./(Dn.*WnS_temp+Ln_temp);
    end
    

if Pb(ii) >= Gar2SpTransition_pressure %Garnet field
    %Fe-Mg Kds  
    if Pb(ii) > 27
    XOlK = 0.35;
    else 
    XOlK = 0.33;
    end   
    XCpxK = 1.04*XOlK;
    XOpxK = 0.93*XOlK;   
    XGarK = 1.4*XOlK;
    XSpK = 1.5*XOlK;
    
    XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMG,YMG,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
    helpDEBUG(ii,:) = [XMGliq, Frac, Wnriiminus1, WnS_temp, XMG,YMG,XOlK,XCpxK,XOpxK,XSpK,XOW];
    
else if Pb(ii) >= Sp2PlagTransition_pressure  %Spinel field 
        
    %Fe-Mg Kds              
    XOlK = 0.33;   
    XCpxK = 0.94*XOlK;
    XOpxK = 0.97*XOlK;   
    XGarK = 1.4*XOlK;
    XSpK = 1.5*XOlK;
    
    
    XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);   
    helpDEBUG(ii,:) = [XMGliq,  Frac, Wnriiminus1, WnS_temp, XMS,YMS,XOlK,XCpxK,XOpxK,XSpK,XOW];
    else % plagioclase field 
        
    %Fe-Mg Kds          
    XOlK = 0.29;   
    XCpxK = 0.94*XOlK;
    XOpxK = 0.97*XOlK;   
    XGarK = 1.4*XOlK;
    XSpK = 1.5*XOlK;
    
    XMGliq = BG15_Mgnum_MassBalance(Frac, Wnriiminus1,WnS_temp, XMP,YMP,XOlK,XCpxK,XOpxK,XSpK,XGarK,XOW);
    helpDEBUG(ii,:) = [XMGliq,  Frac, Wnriiminus1, WnS_temp, XMP,YMP,XOlK,XCpxK,XOpxK,XSpK,XOW];
    end
end



%Plag field
XCP = [XMP(1:4)]*DMP; %bulk partition coefficient in solid
XPP = [YMP(1:4)]*DMP;  %bulk partition coefficient in melt
DP = (XCP - XPP.*Frac)./(1-Frac); 

%Spinel field
XCS = (DMS'*[XMS(1:3),XMS(5)]')'; %bulk partition coefficient in solid
XPS = (DMS'*[YMS(1:3),YMS(5)]')'; %bulk partition coefficient in melt
DS = (XCS - XPS.*Frac)./(1-Frac); 

%Garnet field
XCG = [XMG(1:3),XMG(6)]*DMG; %bulk partition coefficient in solid
XPG = [YMG(1:3),YMG(6)]*DMG; %bulk partition coefficient in melt    
DG = (XCG - XPG.*Frac)./(1-Frac); 


if ii >1
    Wnriiminus1 = Wnr(ii-1);
else
    Wnriiminus1 = Wo;
end

% Calculate concentration of oxide in melt (wt%)
%Plagioclase
XLFP(1) = XOW(2).*Wnriiminus1/(DP(1).*WnS_temp+Ln_temp); %Ti
XLFP(2) = XOW(4).*Wnriiminus1/(DP(2).*WnS_temp+Ln_temp); %Cr
XLFP(3) = XOW(8).*Wnriiminus1/(DP(3).*WnS_temp+Ln_temp); %K
XLFP(4) = XOW(9).*Wnriiminus1/(DP(4).*WnS_temp+Ln_temp); %Na

%Spinel
XLFS(1) = XOW(2).*Wnriiminus1/(DS(1).*WnS_temp+Ln_temp); %Ti
XLFS(2) = XOW(4).*Wnriiminus1/(DS(2).*WnS_temp+Ln_temp); %Cr
XLFS(3) = XOW(8).*Wnriiminus1/(DS(3).*WnS_temp+Ln_temp); %K
XLFS(4) = XOW(9).*Wnriiminus1/(DS(4).*WnS_temp+Ln_temp); %Na

%Garnet
XLFG(1) = XOW(2).*Wnriiminus1/(DG(1).*WnS_temp+Ln_temp); %Ti
XLFG(2) = XOW(4).*Wnriiminus1/(DG(2).*WnS_temp+Ln_temp); %Cr
XLFG(3) = XOW(8).*Wnriiminus1/(DG(3).*WnS_temp+Ln_temp); %K
XLFG(4) = XOW(9).*Wnriiminus1/(DG(4).*WnS_temp+Ln_temp); %Na    



[PCA, DCA] = calculatePCA(Pb(ii), XMGliq, Gar2SpTransition_pressure, Sp2PlagTransition_pressure);
QCA = PCA/DCA;  %Ca content (XLFS)

   
% Calculate NaK# wt, Mg# mol, assign wt% TiO2, K2O, Na/Na+Ca wt
if Pb(ii) >= Gar2SpTransition_pressure
    %Garnet Field
    WNA=(XLFG(3)+XLFG(4))/(XLFG(3)+XLFG(4)+QCA); %NaK# = (Na + K) / (Na + K +Ca)
    WNACA=XLFG(4)/(XLFG(4)+QCA); %Na# = Na / (Na + Ca)
    % assign wt.% TiO2 and K2O
    WT=XLFG(1);
    WK=XLFG(3);
    
else if Pb(ii) >= Sp2PlagTransition_pressure
    %Spinel Field
    WNA=(XLFS(3)+XLFS(4))/(XLFS(3)+XLFS(4)+QCA); %NaK# = (Na + K) / (Na + K +Ca)
    WNACA=XLFS(4)/(XLFS(4)+QCA); %Na# = Na / (Na + Ca)
    % assign wt.% TiO2 and K2O
    WT=XLFS(1);
    WK=XLFS(3);
else
    %Plag Field
    WNA=(XLFP(3)+XLFP(4))/(XLFP(3)+XLFP(4)+QCA);
    WNACA=XLFP(4)/(XLFP(4)+QCA);
    % assign wt.% TiO2 and K2O
    WT=XLFP(1);
    WK=XLFP(3);
end
end

% Calculate mineral components and T (V. 8 MIT+CMASN)
%   Coefficients from Till et al. [2012]


mgratio(ii,:) = [((XOW(6)/40.311)/(XOW(6)/40.311+XOW(5)/71.846)), XMGliq];
NaKratio(ii) = ((XOW(9)+XOW(8))/(XOW(9)+XOW(8)+XOW(7)))./WNA;
%importantComps(ii,:) = [XMGliq XOlK XMGOL (XOW(6)/40.311)/(XOW(6)/40.311+XOW(5)/71.846) (XOW(9)+XOW(8))/(XOW(9)+XOW(8)+XOW(7))./WNA];




water(ii) = CL(H2O_index)*10^-4;



if Pb(ii) >= Gar2SpTransition_pressure
    %Garnet Field    
    
 
% normalized to the 7 tormey mineral components, trying with WNACA
%     Ts(ii)=1302.315+8.768*Pb(ii)-190.805*(1-XMGliq)-477.441*WNA+30.671*WT +469.144*WNACA; 
%     S(1)=0.139+0.006*Pb(ii)+0.254*(1-XMGliq)-0.469*WNA+0.001*WT + 0.346*WNACA ;  %Olivine
%     S(2)=0.177+0.000*Pb(ii)+0.083*(1-XMGliq)-0.365*WNA-0.016*WT + 0.341*WNACA;  %CPX    
%     S(3)=0.693-0.008*Pb(ii)-0.320*(1-XMGliq)+0.338*WNA-0.013*WT+ 0.898*WNACA;  %Plag    
%     S(4)=-0.016+0.002*Pb(ii)-0.011*(1-XMGliq)-0.163*WNA+0.020*WT- 0.318*WNACA; %Qtz   
% %     

%     % normalized to the 7 tormey mineral components
    Ts(ii)=1320.279+8.468*Pb(ii)-203.208*(1-XMGliq)-131.937*WNA+32.949*WT; 
    S(1)=0.152+0.006*Pb(ii)+0.245*(1-XMGliq)-0.214*WNA+0.001*WT;  %Olivine
    S(2)=0.190+0.000*Pb(ii)+0.074*(1-XMGliq)-0.114*WNA-0.014*WT;  %CPX
    S(3)=0.727-0.008*Pb(ii)-0.343*(1-XMGliq)+0.323*WNA-0.008*WT;  %Plag
    S(4)=-0.028+0.002*Pb(ii)-0.002*(1-XMGliq)-0.397*WNA+0.018*WT; %Qtz
         
    
% %    normalized to the 7 tormey mineral components: some reason this
% %    changed from 9/12 to 10/6
%     Ts(ii)=1320.441+8.466*Pb(ii)-203.434*(1-XMGliq)-132.087*WNA+32.954*WT; 
%     S(1)=0.151+0.006*Pb(ii)+0.245*(1-XMGliq)-0.213*WNA+0.001*WT;  %Olivine
%     S(2)=0.184+0.000*Pb(ii)+0.083*(1-XMGliq)-0.108*WNA-0.014*WT;  %CPX
%     S(3)=0.725-0.008*Pb(ii)-0.341*(1-XMGliq)+0.325*WNA-0.009*WT;  %Plag
%     S(4)=-0.026+0.002*Pb(ii)-0.005*(1-XMGliq)-0.399*WNA+0.018*WT; %Qtz
%     
%     %what we were using as of 9/12/17
%     %normalized to the 7 tormey mineral components
%     Ts(ii)=1316.408+8.606*Pb(ii)-199.790*(1-XMGliq)-136.165*WNA+33.326*WT;
% 
%     S(1)=0.167+0.006*Pb(ii)+0.235*(1-XMGliq)-0.215*WNA+0.000*WT;  %Olivine
%     S(2)=0.124+0.002*Pb(ii)+0.125*(1-XMGliq)-0.104*WNA-0.014*WT;  %CPX
%     S(3)=0.735-0.009*Pb(ii)-0.348*(1-XMGliq)+0.329*WNA-0.009*WT;  %Plag
%     S(4)=0.009+0.001*Pb(ii)-0.029*(1-XMGliq)-0.403*WNA+0.018*WT; %Qtz
    
    
    
    Tsdry(ii) = Ts(ii);
    Sdry =S;
    Ts(ii) = Ts(ii) - 26.3619.*water(ii);
    S(1)=S(1)-0.0029.*water(ii);%Olivine
    S(2)=S(2)-0.0093.*water(ii); %CPX
    S(3)=S(3)+0.0075.*water(ii);  %Plag
    S(4)=S(4)+0.0068.*water(ii); %Qtz
    


%     %normalized to the 4 tormey mineral components
%     Ts(ii)=1316.408+8.606*Pb(ii)-199.790*(1-XMGliq)-136.165*WNA+33.326*WT;
%     S(1)=0.147+0.006*Pb(ii)+0.248*(1-XMGliq)-0.087*WNA+0.003*WT;  %Olivine
%     S(2)=0.114+0.002*Pb(ii)+0.131*(1-XMGliq)-0.042*WNA-0.013*WT;  %CPX
%     S(3)=0.727-0.009*Pb(ii)-0.353*(1-XMGliq)+0.603*WNA-0.009*WT;  %Plag
%     S(4)=0.012+0.001*Pb(ii)-0.026*(1-XMGliq)-0.475*WNA+0.020*WT; %Qtz
    
%     %fits using spreadsheet and matlab code (Tim's tormey)
%     Ts(ii)=1316.411+8.605*Pb(ii)-200.055*(1-XMGliq)-135.983*WNA+33.370*WT;
%     S(1)=0.167+0.006*Pb(ii)+0.234*(1-XMGliq)-0.211*WNA+0.001*WT;  %Olivine
%     S(2)=0.122+0.002*Pb(ii)+0.117*(1-XMGliq)-0.103*WNA-0.012*WT;  %CPX
%     S(3)=0.735-0.009*Pb(ii)-0.349*(1-XMGliq)+0.336*WNA-0.009*WT;  %Plag
%     S(4)=0.009+0.001*Pb(ii)-0.026*(1-XMGliq)-0.406*WNA+0.018*WT; %Qtz
    
    
%     Ts(ii)=1303.365+8.750*Pb(ii)-194.648*(1-XMGliq)-5.544*WNA+30.359*WT-23.764*WK;
%     S(1)=0.156+0.006*Pb(ii)+0.239*(1-XMGliq)-0.102*WNA-0.002*WT-0.020*WK;  %Olivine
%     S(2)=0.114+0.002*Pb(ii)+0.121*(1-XMGliq)-0.024*WNA-0.014*WT-0.014*WK;  %CPX
%     S(3)=0.708-0.008*Pb(ii)-0.338*(1-XMGliq)+0.603*WNA-0.015*WT-0.049*WK;  %Plag
%     S(4)=0.016+0.001*Pb(ii)-0.029*(1-XMGliq)-0.481*WNA+ 0.020*WT+0.014*WK; %Qtz

else if Pb(ii) >= Sp2PlagTransition_pressure    %Spinel Field
        
%     %normalized to the 7 tormey mineral components
    Ts(ii)=1212.464+11.966*Pb(ii)-96.899*(1-XMGliq)-89.350*WNA+2.367*WT; %MITCMAS
    S(1)=0.125+0.008*Pb(ii)+0.197*(1-XMGliq)-0.260*WNA-0.023*WT;  %Olivine % All Spinel experiments
    S(2)=0.181+0.000*Pb(ii)+0.071*(1-XMGliq)-0.255*WNA+0.005*WT;  %CPX
    S(3)=0.516-0.001*Pb(ii)-0.130*(1-XMGliq)+0.782*WNA-0.035*WT;  %Plag
    S(4)=0.149-0.005*Pb(ii)-0.098*(1-XMGliq)-0.379*WNA+0.025*WT; %Qtz % All Spinel experiments
    
    
    

    Tsdry(ii) = Ts(ii);
    Sdry =S;
    Ts(ii) = Ts(ii) - 26.19.*water(ii);
    S(1)=S(1)-0.0029.*water(ii);%Olivine
    S(2)=S(2)-0.0093.*water(ii); %CPX
    S(3)=S(3)+0.0075.*water(ii);  %Plag
    S(4)=S(4)+0.0068.*water(ii); %Qtz
%     %normalized to the 4 tormey mineral components
%     Ts(ii)=1212.464+11.966*Pb(ii)-96.899*(1-XMGliq)-89.350*WNA+2.367*WT;
%     S(1)=0.128+0.008*Pb(ii)+0.194*(1-XMGliq)-0.243*WNA-0.017*WT;  %Olivine
%     S(2)=0.186+0.000*Pb(ii)+0.066*(1-XMGliq)-0.242*WNA+0.009*WT;  %CPX
%     S(3)=0.533-0.003*Pb(ii)-0.160*(1-XMGliq)+0.868*WNA-0.015*WT;  %Plag
%     S(4)=0.153-0.005*Pb(ii)-0.100*(1-XMGliq)-0.382*WNA+0.024*WT; %Qtz
      
        
%       %original Till 2012 
%     Ts(ii)=1212.+11.99*Pb(ii)-97.33*(1-XMGliq)-87.76*WNA+3.44*WT-4.58*WK;
%     S(1)=.123+.0085*Pb(ii)+.203*(1-XMGliq)-.268*WNA-.019*WT+.012*WK;  %Olivine
%     S(2)=.180+.0005*Pb(ii)+.075*(1-XMGliq)-.264*WNA+.007*WT+.012*WK;  %CPX
%     S(3)=.554-.0043*Pb(ii)-.192*(1-XMGliq)+.949*WNA-.007*WT-.041*WK;  %Plag
%     S(4)=.143-.0046*Pb(ii)-.085*(1-XMGliq)-.416*WNA +.020*WT+.018*WK; %Qtz
else
    %Plag Field
    
    %normalized to the 7 tormey mineral components
    Ts(ii)=1215.133+10.576*Pb(ii)-72.533*(1-XMGliq)-198.174*WNA+23.590*WT;
    S(1)=0.134+0.006*Pb(ii)+0.171*(1-XMGliq)-0.249*WNA-0.012*WT;  %Olivine
    S(2)=0.230-0.003*Pb(ii)-0.057*(1-XMGliq)-0.139*WNA-0.006*WT;  %CPX
    S(3)=0.390+0.020*Pb(ii)-0.075*(1-XMGliq)+0.485*WNA-0.032*WT;  %Plag
    S(4)=0.225-0.018*Pb(ii)-0.003*(1-XMGliq)-0.229*WNA+0.009*WT; %Qtz
    
        
    Tsdry(ii) = Ts(ii);
    Sdry =S;
    Ts(ii) = Ts(ii) - 26.19.*water(ii);
    S(1)=S(1)-0.0029.*water(ii);%Olivine
    S(2)=S(2)-0.0093.*water(ii); %CPX
    S(3)=S(3)+0.0075.*water(ii);  %Plag
    S(4)=S(4)+0.0068.*water(ii); %Qtz
%     %normalized to the 4 tormey mineral components
%     Ts(ii)=1215.133+10.576*Pb(ii)-72.533*(1-XMGliq)-198.174*WNA+23.590*WT;
%     S(1)=0.135+0.006*Pb(ii)+0.166*(1-XMGliq)-0.240*WNA-0.004*WT;  %Olivine
%     S(2)=0.234-0.004*Pb(ii)-0.056*(1-XMGliq)-0.126*WNA-0.002*WT;  %CPX
%     S(3)=0.398+0.017*Pb(ii)-0.116*(1-XMGliq)+0.567*WNA+0.001*WT;  %Plag
%     S(4)=0.233-0.019*Pb(ii)+0.007*(1-XMGliq)-0.202*WNA+0.005*WT; %Qtz
    
    
    %original Till 2012
%     
%     Ts(ii)=1216.+10.44*Pb(ii)-72.83*(1-XMGliq)-194.9*WNA+24.08*WT-1.55*WK;
%     S(1)=.132+.0069*Pb(ii)+.173*(1-XMGliq)-.261*WNA-.009*WT+.010*WK;  %Olivine
%     S(2)=.238-.0051*Pb(ii)-.065*(1-XMGliq)-.097*WNA+.004*WT-.013*WK;  %CPX
%     S(3)=.408+.0146*Pb(ii)-.136*(1-XMGliq)+.630*WNA+.014*WT-.028*WK;  %Plag
%     S(4)=.222-.0164*Pb(ii)-.029*(1-XMGliq)-.272*WNA -.009*WT+.031*WK; %Qtz
    end
end

testPboundary_Sp2Gar = (1.2  - .52*WNA+2.11*XMGliq)*10;
testPboundary_Plag2Sp = (0.823  +1.138*WNA+0.165*XMGliq)*10;

S = S./sum(S);
Sdry = Sdry./sum(Sdry);
%test(ii,:) = [Pb(ii) Ts(ii) Tb(ii) (1-XMGliq) WNA WT WK testPboundary_Plag2Sp testPboundary_Sp2Gar QCA];
test(ii,:) = [Pb(ii) Ts(ii) Tb(ii) (1-XMGliq) WNA WT WK WNACA testPboundary_Plag2Sp testPboundary_Sp2Gar QCA 1 1 ];
testS(ii,:) = [S sum(S)];
testSdry(ii,:) = Sdry;



end




